

















Vo.tumE 18 MARCH, 1951 NuMBER 3 





CONTENTS 


Turbulent Boundary Layer in Compressible Fluids E. R. VAN DRIEST 


First- and Second-Order Theory of Supersonic Flow Past Bodies of Revolution... 
Of SAREE he SAE a Pes a ae A ee eM er. aes? MILTON D. VAN DYKE 


Methods for Calculating the Flow in the Trefftz-Plane Behind Supersonic Wings 
P. A. LAGERSTROM AND MARTHA E. GRAHAM 


Experimental and Theoretical Studies of Flame-Front Stability 
GEorRGE H. MARKSTEIN 


Readers’ Forum: ‘Thickness of a Steady Shock Wave,” by J. J. Bernard; ‘‘Com- 
ments on ‘A Type of Lifting Rotor with Inherent Stability,’” by G. J. Sissingh; 


“Buckling of a Pretwisted Pin-Ended Column,” by G. V. R. Rao; ‘On the 
Equations of Longitudinal Stability—II,” by John W. Miles; ‘Note on the 
Relation of Lifting-Line Theory to Lifting-Surface Theory,” by Eric Reissner; 
“The Application of Nonstationary Airfoil Theory to Dynamic Stability 
Calculations,” by E. V. Laitone and E. R. Walters 210-216 


Published by the 
INSTITUTE OF THE AERONAUTICAL SCIENCES, INC. 











Associate Editor: 


Aerodynamics 


Francis Clauser 
Howard Emmons 
R. Paul Harrington 
R. T. Jones 

Carl Kaplan 

W. B. Klemperer 
John G. Lee 

H, W. Liepmann 
Cc. C. Lin 

H. F. Ludloff 
John R, Markham 
Clark B. Millikan 
Shatswell Ober 

W. Bailey Oswald 
Allen E. Puckett 
Elliott G. Reid 

H. J. F. Reid 
Russell Robinson 
G. B. Schubauer 
W. R. Sears 

John Stack 
Theodore von K4rman 
Richard von Mises 


Flight Testing 


W. K. Ebel 

R. R. Gilruth 
Melvin N. Gough 
Robert Stanley 
E. G. Stout 


Meteorology 


C. F. Brooks 
D. M. Little 
E. J. Minser 
F. W. Reichelderfer 
H. C. Willett 


JOURNAL OF THE 
AERONAUTICAL SCIENCES 





Editor: Hucu L. Drypen, Director, N.A.C.A. 
Managing Editor: Berneice H. Jarcx 


Rosert R. DEXTER 





Editorial Committee 


Structures and Materials 


W. B. Bleakney 
E. W. Conlon 

J. P. Den Hartog 
L. H. Donnell 
E. C. Hartmann 
. J. Hoff 

. E. Lundquist 
. F. Nagel, Jr. 
8 Newell 
fred S. Niles 

. Ramberg 

. Reissner 

. R. Shanley 

1. W. Sibert 

C. R. Strang 

S. Timoshenko 
H. S. Tsien 
Chi-Teh Wang 


Rape Cons 


Rotating Wing Aircraft 


Robert Coleman 
K. H. Hohenemser 
Bartram Kelley 
Ralph McClarren 
R. H. Miller 

R. H. Prewitt 
Igor I. Sikorsky 


Physiologic Problems 


H. G. Armstrong 

Louis H. Bauer 

Eugene Du Bois 

W. Randolph Lovelace, II 
Ross A. MacFarland 

A. D. Tuttle 


Flight Propulsion 


George W. Brady 


Frank W. Caldwell 


L. S. Hobbs 

E. E. Johnson 
Joseph Keenan 
R. P. Kroon 
Arthur Nutt 
Hans Reissner 
Abe Silverstein 
C. Fayette Taylor 
E. S. Taylor 

P. Taylor 

E. S. Thompson 
H. S. Tsien 


Aircraft Design 


W. E. Beall 

W. K. Ebel 

Hall L. Hibbard 
Richard Hutton 
C. L. Johnson 

I. M. Laddon 

C. J. McCarthy 
Fred E. Weick 
Robert J. Woods 


Air Transport 


Charles Froesch 
J. A. Herlihy 

R. D. Kelly 
Otto E. Kirchner 
Jerome Lederer 
John C. Leslie 
R. Dixon Speas 





Instruments 


Allan G. Binnie 
W. G. Brombacher 
Charles H. Colvin 
B. M. Craig 

B. Del Mar 

C. S. Draper 

O. E. Esval 

W. L. Howland 

J. E. Lindberg 
David W. Moore 
George H. Purcell 
W. A. Reichel 

C, F. Savage 

O. Hugo Schuck 
C. L. Seward, Jr. 
R. Sylvander 

D. K. Warner 

P. F. Weber 


Vibration and Flutter 


Lee Arnold 

M. A. Biot 

William Bollay 
Chieh-Chien Chang 
Martin Goland 

J. P. Den Hartog 


Fuels and Oils 


D. P. Barnard 
J. H. Doolittle 
Graham Edgar 
R. T. Goodwin 
S$. D. Heron 

A. M. Rothrock 





SUBSCRIPTION RATES 


Journat or THE AzRonauticat Sciences, Published Monthly.—United States and Possessions: 1 Year, $12.00; Single Copies, 
$1.50. Foreign Countries Including Canada (American Currency Rate): 1 Year, $13.00; Single Copies, $1.50. 

Aeronautica. Enorneertne Review, Published Monthly.—United States and Possessions: 1 Year, $3.00; Single Copies, $0.50. 
Foreign Countries Including Canada (American Currency Rates): 1 Year, $3.50; Single Copies, $0. 50. 

Notices of change of address should be sent to the Institute at least 30 days prior to actual change of address. 

Manuscripts for publication, proofs, and all correspondence should be addressed to the Editorial Office of the Journat. 
Correspondence regarding subscriptions may be addressed to Publication Office, 20th and Northampton Sts., Easton, Pa., or to the 
Editorial Office of the Journat or THE AERONAUTICAL Sciences, 2 East 64th Street, New York 21, New York. 


Copyright, 1950, by the Institute of the Aeronautical Sciences, Ine. 


Entered as Second Class Matter at the Post Office, Easton, Pa., May 1, 1937. Acceptance for mailing at a special 
rate of Postage provided for in the Act of August 24, 1912. Authorized April 29, 1937. 














JOURNAL OF THE 
AERONAUTICAL SCIENCES 





VOLUME 18 


MARCH, 1951 


NUMBER 3 





Turbulent Boundary Layer 
Compressible Fluids 


In 


E. R. VAN DRIEST* 


North American Aviation, Ine. 


SUMMARY 


The continuity, momentum, and energy differential equations 
for turbulent flow of a compressible fluid are derived, and the 
apparent turbulent stresses and dissipation function are identi- 
fied. A general formula for skin friction, including heat trans- 
fer to a flat plate, is developed for a thin turbulent boundary 
layer in compressible fluids with zero pressure gradient. Curves 
are presented giving skin-friction coefficients and heat-transfer 
coefficients for air for various wall-to-free-stream temperature 
ratios and free-stream Mach Numbers. 

In the special case when the boundary layer is insulated, this 
general formula yields skin-friction coefficients higher than those 
given by the von Karman wall-property compressible-fluid for- 
mula but lower than those given by the von Karman incom- 
pressible-fluid formula. Heat transfer from the boundary layer 
to the plate generally increases the friction and heat-transfer 
coefficients. 


NOMENCLATURE 


Variables 
time 
coordinate along plate in direction of free stream, 
measured from the leading edge 
coordinate normal to plate, measured from the 
plate 
2 = coordinate along plate norma! to free stream 
u,v x and y velocity components, respectively 
Pr, Py, pz normal stresses in the direction of x-, y-, and z- 
axes, respectively 
shear stress 
shear stress in x-direction in plane normal to 
y-axis 
shear stress in y-direction in plane normal to 
x-axis 
heat transfer 
heat transfer in x-direction 
heat transfer in y-direction 
fluid density 
coefficient of viscosity 
coefficient of heat conductivity 


Presented at the 1950 Heat Transfer and Fluid Mechanics 
Institute, Los Angeles, June 28-30, 1950. Received April 28, 
1950. 
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absolute temperature 
specific heat at constant pressure 
eddy coefficient of viscosity 
eddy coefficient of heat conductivity 
ratio of specific heats 
Mach Number 
mixing length 
Sutherland constant 
thickness of the boundary layer when the fluid 
is compressible 
thickness of the boundary layer when the fluid is 
considered incompressible 
Cy = total coefficient of friction 
aa local coefficient of friction 
R Reynolds Number based on length x 
Npr = Prandtl Number 


Su scripts 

w = wall condition 

x free-stream condition 
laminar 
turbulent 


lam. 
turb. 


INTRODUCTION 


aoe MAJOR PROBLEMS ENCOUNTERED TODAY in 
aeronautics are the determination of skin friction 
and skin temperatures of high-speed aircraft. Since 
the friction drag is a considerable portion of the total 
drag of a guided missile, it follows that miscalculation 
of the friction drag can result in considerable error in 
missile range. Furthermore, skin temperature is a 
decisive factor in the structural design of a high-speed 
missile. 

These two problems arise as a result of the presence 
of the fluid boundary layer. Whether the boundary 
layer on a given missile is laminar or turbulent under 
certain conditions is as yet uncertain. However, in 
the present paper only the turbulent case is discussed. 
In particular, the purpose of the paper is to derive a 
general formula for skin friction including heat transfer 
to a flat plate for a thin fully turbulent boundary layer 
in compressible fluids with zero pressure gradient 
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This work is a combination of two previous reports! * 
by the same author. 


A first attempt to obtain theoretical results on tur- 
bulent skin friction on an insulated flat plate was made 
by von Karman when he suggested that the fluid state 
at the wall be used in the incompressible fluid formula; 
it seems that this suggestion is too optimistic because of 
the fact that there is a strong variation in temperature 
throughout the boundary layer. Frankl and Voishel‘ 
used von Karman’s similarity law for mixing length and, 
as a result of compressibility effects, encountered 
considerable complication in deriving an expression for 
turbulent friction with or without heat transfer. 
Ferrari’ has made a general study of the problem; 
however, his apparent turbulent shear stress for com- 
pressible fluids is not in agreement with that of the 
author, and no engineering formula for skin friction 
or heat transfer is presented. Wilson'® recently de- 
rived a formula for insulated plates using the von 
Karman similarity law; his formula is similar to, but 
not the same as, the formula given in this paper. 


In the present paper the differential equations for 
turbulent flow of a compressible fluid are derived for 
the purpose of identifying the apparent turbulent 
stresses and dissipation function. The boundary layer 
is then assumed thin, and the usual relation [Eq. (34) ] 
between temperature and velocity is obtained when the 
turbulence Prandtl Number ce/x is taken as unity, 
where c, is the specific heat at constant pressure, ¢ is 
the coefficient of eddy viscosity, and x is the coefficient 
of eddy heat conductivity. The Prandtl wall differ- 
ential equation is next derived for compressible fluids— 


di _ bite 
dy KYpy 


Here, u% is the mean local velocity, 7, is the mean shear 
at the wall, j is the mean local fluid density, y is the 
distance from the wall, and K is the proportionality con- 
stant in the mixing length formula/ = Ky. Finally, 
a general formula for skin friction including heat 
transfer is developed, and curves are presented giving 
skin-friction coefficients and heat-transfer coefficients 
for air for various wall-to-free-stream temperature 
ratios and free-stream Mach Numbers. For the in- 
sulated plate case, it is seen that the writer’s results 
for skin friction lie above those of von Karman. Cool- 
ing the wall generally increases the skin-friction and 
heat-transfer coefficients. 


Viz., 


THE THIN TURBULENT BOUNDARY LAYER 


The basic differential equations describing the two- 
dimensional flow in a boundary layer on a flat plate are 


ra) 
Continuity: eS + — > (ou) +5 >. (ot) = (1) 
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po Op: . d 
“ns pus — - + os oe + — es 
oy Ox oy 
Momentum: > dp > 
4 a ee ae ee 
Pat - mus Ox * ¥ Oy oy - Ox (2b) 
| oO oO 
\p - (¢pT) + pu = > (GD) + oka (GT) - 
a a A 
Ot Ox eg Ox - 
Energy :1 : 
(p + Peo" + (+ p)2 > + 
Ou Ov 
( ” oy lias Ox 9) 
where 
p = —(1/3) (be + py + p:) (4) 
4 —F a (Se + oe) + 9, Ou P 
Pra itis Ox °) 
4 2 » (OX + ae s+ ow 
b+ by = ~3H + 2u > (6) 
i _2 (= =) (7 
Ore ore ee ‘ 
+2) 
7” —e Oy Ox °) 
oT oT 
1 = kT, k— 
$ >" oy 


The last four terms on the right-hand side of Eq. (3) 
are called the dissipation function because they are 
functions of the viscosity. 


The Mean Momentum Equations—Reynolds Stresses 


The above momentum equations are rewritten in the 
following form by using the continuity equation: 


ve (pu) +2 (pu?) +2 _—— = 2 es (9a) 
re) OTzy 

> (oe) += (pu) += > (on aay Se 
oy ox 


The following substitutions are then introduced into 
Eq. (9a): 


Tyz + Tos 


+ pr 


u=utu’, p=pto, m= 


pu = pu + (pu)’, pv = pu + (pv)’, pr = pz 


where the bars indicate slowly varying temporal meat 
values and the primes indicate instantaneous fluctua- 
tions from the mean. Time average of the resulting 
equation yields 
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J. (3) 
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ing 





TURBULENT BOUNDARY 


°F ii + p’u’) + 5 Lon: u+ (pu)’u | + 
0 a ae = Op: OT yz 
2 iat Gaal = 242 9 


By time average is meant, for example, 


7 2) 
i= — u dt 
TY). - (7 2) 


also, it can be shown that 0/0 = 0-/2. 


The equation of continuity, when averaged, becomes 


+> Glt= @< 11) 
aT he (pu ay (pv) = ( 


With some rearrangement and use of Eq. (11), Eq. (10) 


then reduces to 


. 


LAYER IN COMPRESSIBLE 


FLUIDS 147 


D re) 
+ pv (p’v’) + 


or dx dy ar 
o = ——— 
[Py — (pv)’v pv) ’v’| *; a — (pu)'v’| (12b) 
oy 


The boundary layer is now assumed to be thin— 
i.e., 7 << a. Hence, it is seen that terms on the 
left-hand side of Eq. (12b) are negligible compared to 
the corresponding terms in Eq. (12a), so that Eq. 
(12b) reduces to an equation devoid of mean motion 
terms. Now, in Eq. (12a), 0p,/dx is to be dropped 
because any external pressure gradient is omitted. 
(0/dx) [—(pu)’u’] should be negligible compared to 
(0/dy) [—(pv)’u’| because of thinness of the layer. 
Finally, 07,,/0y definitely can be dropped because Ty; 
is the average of the viscous shear, analogous to in- 
compressible fluid experience. The resulting momen- 
tum equation involving mean motion terms for thin 
turbulent boundary layers is, then, 














Pa we Ou e+ —Ou% oO (—yul) + 
= pares - 
ro oe Oy oe te 
- y = — _— a] 
Qo — a? . Ot = ox - O% ot ane oy dae 
oe [p: — (pu)’u’| vS > Fr — (pv)'u’] (12a) (13) 
Comparison of Eqs. (12a) and (2a) shows that if, in and in the mean steady state 
Eq. (2a), p, u, pu, pv, pz, and 7,2 are replaced by their da da _ 
mean values, then the terms (0/Ot) (—p’u’), (0/Ox) X pu “a pu rm oo ov [—(pv)’u’] (14) 
{—(pu)’u’|, and (0/dy) [—(pv)’u’] must be added to the , ‘ 
right-hand side of the equation. The terms —(pu)’u’ Nowlet 
eae ’ > > > ‘ ¢ . 
and —(pv)’u’* will be called the Reynolds or apparent ~igu)’u! = €(0a/2y) (15) 


pressure and shear stresses, respectively, in the case 











of the turbulent flow of a compressible fluid. In the where e« shall be called the eddy viscosity as in incom- 
same manner, Eq. (9b) reduces to pressible fluid theory. Eq. (14) then becomes 
Se “—s aia: lia tmx <a Sake ws Ou Ou re) ( 2) 
_* Expanded, this term becomes (pv)’u’ = pu'v a 16 
dp'u’ — p'u'v’. ” Ox + oy oy : ov ats 
The Mean Energy Equation—Apparent Dissipation Function 
Use of a (1) and (2) allows Eq. (3) to be written in the form 
) u re) us re) puu? re) 0 /puv? 
(pcpT) oe S. (pucsT) +5 > (p0esT) +3 os “) +— - + (' - + — (= + — ( : + 
ar Ox oy 2 ot ox \ 2 
re) “' Op Op pe Ody +4 Op: ore Op, oy 
at Sew Fane ae Se oe oe } = — - + -~** + ys. t 
oy ( ot . ox oy Ox ” oy lies oy st Ox aii » + bs) x 
Ov Law 
(p + py) .* re a+ "es (17) 
As before, the following substitutions are then introduced: 
u=atu T T+ T’ p =pt+p’ 
v mg: p =ptp’ be = pr + Pr’ 
pu = pu Sua Tye = Tye + Tyr = Py = Py + Py’ 
pu = pu + (pv)’ Try = Ty + Try" Ww = qv 7 qv’ 
Qe = G+ a’ 


Where the mean values again are slowly varying. With these substitutions inserted, Eq. (17) becomes 








ra) e- 
© lal +O + TH+; bey Low + (pu)'| (T+ T)}+ > feoloo + (o0)'] (T + T')} + 


: K (p+ p’) (a+ “ae of ef lou + (pu)’| (a + ure + de b [pv + (pv)’| (a + u'y2 + 
or 7 3 p p ( dy pe pu f 
ra) l 2 fe) , x 0 oat? a Prd ol _ 
5 [5 e+e) @ + ee 15 lou + + (pu)'| (6 + +e; + 15 loo + + (p2)'] (@ + 2’)? = 


fe) / 
(qQy + Q@) + 


fe) / — / fe) , = / ) 
soeneusngornsernzoenezarees 


oO a) 
(i + u’) 9 (Tye + Ter’) + (4 + wy (pz + pr’) + (+ 0’) a (Try + Try’) + (8+ 0 2 (Pu + p,’)+ 


fa) 2 0 2 
(b+ p' + pe + Pr’) Gaitu)+(pt+p’+p,t+ by’). (6 + 0’) + (tyr + Tyr’) Sy (a4 + u’) + 
; fa) ss , O - , o f = , / - \l 
(try + try’) — O@+ 0!) =— (PtP) + — b+ 0’) + (p. + br’) (4H + u’)f + 
Ox ot Ox 


Oo 
(qv + q’) + 5 U(r + tn’) i + 0’) } + 


fe) z , , = nl fe) , 
(b+ b’) (by + by @+0')f +— (a tau’) + 
Oy Oy 


Ox 


Gata OFe) 


Upon multiplication of the terms in the brackets and averaging term by term over time, there results 


o 
> (oT + G p'T’) + 5 © Ley pus T + cy (ou)’T” |+< 1 lee po T + Cy (pv)'T’| + 
lo -~% 9 = / / / 49% = ® fe - / / , 9? 
ome + re + ve + p’u yp Lowi? + puru’® + 2a (pu) + (pu)'u’?| + 
lo -9 12 = ee, , $a,82 
va (pvt? + pu-u’? + 20 (puv)’u’ + (pu)’u’?| + ay (ee D4 By’? + D-p'y’ + p’v’”?| + 
] 3) =9 19 = or, 772 I re) -9 72 s a: 779 
|pu-d? + purv’? + 28-(pu)'v’ + (pu)'v’?| + [pv 3? + pu-v’? + 28 (pv)'v’ + (pv)'v’”?| = 
2 Ox Me 
Op er, rs Og: Oq, 
oP + Sb + pevit + plu’ + piu’) + > (p-d + pd + p'v’ THVT e Te, * 


fe) - , / ) = Fi¥ 
a (tye° + Ty, ') + -. tig Bl A tay V') 


Carrying out the indicated differentiation of Eq. (19), one obtains 


(pu) + pu 


Ox 


_ Op om fe) a) 
Ss + re (cyl ) t 5 (c, p' 7") + ef? “a (cp7 ) +5 > fan ( (pu)’/T’| + qf - aa (pv) + 





pv . (c)T) + 5 ° Io (pv)'T’| + 3 a + . 2 +p a + — + u - (p'u’) + WS + 
. €e * + i + = = ( (pu) + pu oc ‘) +" + . a ea + i 5 =  [ou) ‘0’ + (ou) + 
= Da ls + %- pu "> + * (pv) + po, “oy (3 od dy )+a, > (e0) u’| + (pv) "s + 
: ~ (pu) + pus > yee - ian +70 ° [(pu)’v’| + oes "a 4. 3 (pv) = + 
| (pv) + po 5 ae “ 2a (pv) + 0 : [(ov)’v"] + oy + rs ioe ‘| = oP + i + 
12h 5 + u oP: 4° - (pu) + © Dn (pz w) + BS +0 + Dig tO + OF) + 

: (py’v’) + a + ® + a + — + 5 (ty2'u’) + rao + 6 as + ~ (T2y'v" 
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(18) 


(19) 


(20) 
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This equation can be reduced by use of the continuity and momentum equations. The continuity equation |Eq. 
(11)] eliminates the Ist, 4th, 7th, 11th, 13th, 18th, 20th, 25th, 27th, 32nd, 34th, 39th, 41st, 46th, and 48th terms 
on the left-hand side. The momentum equations [Eqs. (12a) and (12b)] eliminate the 10th, 14th, 17th, 21st, 24th, 
98th, 31st, 35th, 38th, 42nd, 45th, and 49th terms on the left-hand side and the 5th, 11th, 17th, and 20th terms on 
the right-hand side of the equation. Collecting the remaining terms then gives 

7, (¢p7) + pu 2° (cpT) + pv . , (6 yl) — oe ie i oF d ew ~ (~¢,9'T") — pu -” ia - 


ot Ox Oy 
ei? fe” | = Oo (u Se = a = 
BE) ast) a8 (0 2) a (0 + re) + 


-+4 
0 ) | owe -) 2 [ey v) 2 a") 
oe x ) a ly! /'v' — * ly! a 
fa (p'u' + prt “a x + a (p’v’ + p,'v’) a “(7 ) ze ; 4 


oO — O | (pv)’u’? Og: dy t 
‘ (tyz't’) — 2| 9 | + dx > [es(ou) ‘F*) +° oa is > ti Cp(pv)’T’ | + (pb + py) — 
vo ee oe (ao? , OV 4 dit Tag , On 4 o# 7 , 00 (21) 
1 — (pv)'v wi (ee  * v 2 
wie Pol ay — SVT aT te a ee OT Pg MY oe 


Comparison of Eqs. (21) and (3) shows what terms must be added when the instantaneous values of the variables 
are to be replaced by their mean values. The last eight terms of Eq. (21) correspond to the last four of Eq. (3); 
hence, the additive terms 

—(pu)’u’ (0% /Ox), —(pv)'v’ (08/Ov), —(pv)’u’ (O%/OyY), —(pu)’v’ (00/Ox) 
will be called apparent dissipation function terms. 

Eq. (21) can be reduced further. (1) In Eq. (19), wu’? and vy’? are expected to be small compared to “°, and, since 
iis assumed of the same order as c,7, they can be dropped from the equation. (2) The triple correlations are also 
expected to be small and, hence, can be neglected. Next, because of thinness of boundary layer and absence of 
external pressure gradient, terms containing pressure gradients should be negligible. Hence, Eq. (21) becomes 


(cpT) + pu ~ (cp) + pv = (cy) _ oe = ~ (—c, p'T’) — pn’ = — p’v’ 4 + 
° a) +2 Gie) + +2 |-. YP) + 2 + 2 [-c Tl + GO + pa) - 
aid wy” ox OX oy oy ° 7 
— Of : ey — Ow On _— On Ov 7% oy 
(pu)'u = + (p+ p,) i — (pv)'v ‘* + Tr oy — (pv)'u re + Tey dy — (pu)’v’-s— (22 


Terms representing molecular action can further be and in the mean steady state 
neglected compared to terms representing molar action; > 


: ; a) o om ; 
hence, the 4th, 5th, 6th, Sth, 14th, and 16th terms on — py (cpl) + pv (cpl) = | —cp(pv)'7”| ae 
the right-hand side can be ignored. Because of thin- Ox oy oy 
ness of the layer, (0/Ox) |—c,(pu)’T’] should be small (pv) 1’ Ou (24) 
; oy 


relative to (0/Oy) |—cp(pv)’T’| and can be dropped. 
Furthermore, 07 /Ox, 00/Oy, and 02/Ox are much smaller 
than 07 /Oy, and, hence, all remaining terms containing 
these as factors will be assumed small and will be ig- 
nored. Also, —p’v’ (O@/Of) is dropped, since it is ex- , . os 

: hint i PI ‘pte where « will be called the eddy heat conductivity as 
pected to be small relative to —p’u’ (07/01) because .. : : : a . 
= ‘ ‘ ; : in incompressible fluid theory. Hence, with Eqs. (15) 
?< a4. As a final approximate result for thin turbu- fee rer : 
le : and (25), Eq. (24) becomes 
ent boundary layers, Eq. (20) reduces to 


The following definition is now introduced— viz. : 


—c,(pv)'T’ = «(OT /dy) (25) 


2 > > : fe) (,T) + ) (c,T) re) ( ~) 4 (o") 
7 7 -" Cp v Cpl) = iN € 
el) +n kl+e— eh -= « oot ee od "ts oy \“ oy oy 
of Ox oy Ot oF, 
ou Oo PY aa 
ay Ot t ol i aladll al oy as ool po)" r | 7 A Solution of the Energy Equation 
(pv)’u’ Ou (23) The mean steady turbulent-flow equations for con- 


Oy tinuity, momentum, and energy for thin boundary 








150 JOURNAL OF THE 


layers are collected 


(0/Ox) (pu) + (0/Oy) (pv) = 0 (27a) 
Ou Ou re) Ou 

, v= — 27b 
ait "hha = =) aia 

ra) ia Ze = ra) ol ou\? 

ae Oe, aS ye 

sl ox ot) te oy (Gt) = > (« af) bi (5*) 

(27c) 


Now, analogous to the procedure with laminar bound- 
ary layers, where the Prandtl Number c,u/k is taken as 
unity because its actual value for air is near unity, a 
turbulence Prandtl Number will be invented which will 
be taken equal to unity. Thus, it is assumed that 
c,e/k = ] (28) 


p*/ 


This assumption seems justified by the fact that « and 
«x result from the same internal mechanism in turbulent 
flow just as uw and k result from the same internal mech- 
anism in laminar gaseous flow. Hence, since c, is 
constant, Eq. (27c) becomes 


(cT) - 


2 on\? 
dy Ee (¢ ws r) | + (5) (29) 


The temperature is now assumed to be a function of 


ra) a re) 
amy cs oe 
us. (cpl) + pr a 


the x-component of velocity—i.e., c,>T = f(a). There- 
fore, Eq. (29) becomes 
— On df(ii) - Ott df ( (a4) _ oO ol. Oi we) 
pu i + 
Ox du ” oy dit Oy du 
sy _ df(i) > (: 2) « (8) SLO 4 
“\ay/ da dy Voy dy) di? 


du\? — df(u%) O (. =) Ee @) 
« (5) di dy\ dy lia di* il oy 


(30) 
Hence, the momentum Eq. (27b) is satisfied if 
d*f(u)/du? = —1 (31) 
which upon integration gives 
Cpl = A; + By — (a?/2) (32) 


where A, and PB, are constants of integration. 

The Prandtl Number of the laminar sublayer is 
next assumed equal to unity. Since in the case of 
laminar flow an equation identical to Eq. (32) is ob- 
tained for Prandtl Number equal to unity,° it follows 
that Eq. (32) may be used across the laminar sublayer 


to the wall. Hence, with the boundary conditions 
T =T,,atuz =0 
IT =T, atu =u, 


and the fact that 


Uo? /Wypl @ = [(y — 1)/2]M..? 
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Eq. (32) becomes 
= —_ Tw Ga m un /- 
i a * ii. 
5 eee ~ u 
— Ng (1 — ) (34) 
2 U co Ue. 


where the subscripts © and w refer to the free stream 
and wall, respectively; y is the ratio of the specific 
heats; and M, is the ratio of the free-stream velocity 
to the velocity of sound in the free stream. 
has previously obtained this equation.) 
Differentiation of Eq. (34) yields 


1 /oT 1 r=) y¥-1 On 
eo -} = — 1— M.? - 
(ge (5) t co  ( ce i 2 is 
it dil 


Uo" OY 


(Croceo! 


(y — 1)M..? 
which becomes, at the wall, 
oT ie ,. —1 i 
( ) ™ (1-72) +7 uw: |(2) 
OV / » Uo t cs 2 Oy), 
(35) 


or, introducing the coefficient of heat conductivity at 
the wall, k,, and the coefficient of viscosity at the wall, 


Muy 
(2 
Des 


oT ae ae 
Rw ( ) = . 
OV) Kw Ua 
y-—1 | 
Bn? th es 
2 |. (3 


whence, by definition of g and 7 and the assumption 


that Cohw, Rw = L 
M . | (36) 


Vw Fe ( z* ) ; aw l 
= ty 1-=")+ 
Tw ” Ue ‘a 2 
valid for turbulent, as well as laminar, flow. 


In the special case when the boundary layer is in- 
sulated, g,. = 0 and Eq. (36) yields 


T/T. = 1+ [(y — 1)/2]M.? (37) 
which, when substituted into Eq. (34), gives 

= = = - a ' Me? () (38) 
Use of Eq. (33) then gives 

T/T. = (Tw/T a) — (#?/2pT w) 
or 

Cpl» = Gl + (2/2) = constant (39) 


Eq. (39) states that, in the case of an insulated flat 
plate, the energy content per unit mass is constant 
across the boundary layer whether laminar or turbulent. 

In the general case of heat transfer to or from the 
boundary layer, Eq. (36) can now be rewritten, with the 
aid of Eq. (37), in the form 
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TURBULENT BOUNDARY 
Qu = (CpTw/Uo) (Twin, — Tw) (40) 
where 7'»,,, is the temperature the wall would acquire 
if the boundary layer were insulated. 

The film coefficient of heat transfer, h, will be defined 
by 


Go = h(Tw,.. — Tw) (41) 
so that 

h = Cp (Tu/Uow) (42) 
Defining a dimensionless heat-transfer coefficient by 

Cx, = h/Cppotta (43) 
and the local coefficient of friction by 

Cr, = Be/Palta® (44) 
it follows that 

Cu. = (1/2)6r (45) 


Therefore, once the local coefficient of friction is known, 
the local coefficient of heat transfer is directly obtained. 
Any refinement in the above procedure would require 
further knowledge of the turbulence mechanism and 
the extent of the laminar sublayer. 


TURBULENT SKIN FRICTION WITH HEAT TRANSFER 


It has been shown in a previous section that the tur- 
bulent shear stress is given by 
r = —(pv)'u' 
ae ee ar ee ae | 
=-—-fuv -ten' ~ ga 
In the case of the thin boundary layer where v ~ 0, 
this equation yields 


r= — pu'v'’ — p'u'd’ (46) 
Now, if the triple correlation is neglected, Eq. (46) re- 
duces further to 


r= —pu'v’ (47) 


which has the same form as for incompressible fluids. 
Introduction of the Prandtl mixing length, /, then 
gives 

rt = pl*(du/dy)? 
Following Prandtl’s method, a wall formula is next ob- 
tained by assuming that / = Ay and that the shear 
stress is constant near the wall; thus, 


Ty = pK*y?(dy/du)? 
or, rearranging, 
dii/dy = (1/K) Vr./p (1/y) (48) 


which is the same as the Prandtl incompressible fluid 
wall formula, except that the density is variable. 

Since the pressure is constant through the thin 
boundary layer, there results from the perfect-gas law 
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v/F 


Now, rearrangement of Eq. (34) yields 


i y—1 \E | u 
= ] 1 — M.? — 1 = 
T. 7 ( ¥ 2 Te U ew 


y-1 et y a 
M.* 50 
2 ; a 5 50) 


B/ Pw = 


so that from Eq. (49) 


Dp l 
- = . — (51) 
Pw 14+ B(a#/u.) — A*(a/u.)? 
where 
— |] — | 
1—— Me? 1+ 7— mM.’ 
4? = ~ and B = - — | 
: — " T./T. 


Substitution of Eqs. (51) and (48) gives 


d(u Uc) _ _ 1 ae dy 
E +B(") -a(2)]" Ua KV pw ¥ 
1 


U co bx 


(52) 


integration of which yields 


; , 2A? (u4/u.) — B 
sin-! —- 


A (B? + 44%)” = constant + 


3 y" (53) 
” a ny Je 
ite" 


Since Eq. (53) must reduce to the incompressible fluid 
case when M,, = Oand 7,,/7.. = 1, it will be necessary 
to write it in the form 


, 2A u.) — B 1 B 


Nalin sin-! —; = 
(B? + 442)” a sin (B? + 443) 


l iu 
constant + -—a— ing 
Uno KY po 


sin 


or 
, 2A? (hi uo.) — B l B 


. 4 1 7 _ 
(B? rn 4A2)' P i A sin (B? + 4A2)'/2 


1 \" (F + 1 1 y" ».) 54) 
a4 n _ e 
Uw Pw K Pw Yw ' 


where F is a constant and v,, the kinematic viscosity at 
the wall, is introduced because of its influence in the 


sin 


laminar sublayer 
Now, the wall shear stress is computed from the for- 


d 5 
Tf. = pu(u. — u)dy (55) 
dx 


in which the double and triple correlations have been 
neglected and where 6 is the thickness of the boundary 
layer. When Eq. (51) and dy obtained from Eq. (54) 
are substituted into Eq. (55), there results 


mula 
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Du,u. d § las : oo B } 
Ty. = —— -—— | Ja? Gxp sin~ - or oo 
SK “de VP LA (BEF 444] 
(56) 
where 
D=e '* 
" / 
a = Ku./V ru/ po 
, | ’ 2(1 — 2) ye 
. o (1 + Bz — A?’s?)'” 
exp k nt a8 — . | Iz (57) 
>, sin az o4 
A (B? + 44%)? 
andz = /u.,. In these equations, x has itsusual mean- 


ing of distance along the flat plate in the direction of 

the free stream and measured from the leading edge. 
Eq. (57) for J can be expanded in a series with the 

aid of integration by parts. The series can be approxi- 


mated by 


l 
= 70 +B A” * 
a [a a 2A? — B l 
sp tt gn | ta 
Ge. B wx 
exp | — r sill (B24 1A2)' ] (58) 


when terms of higher order than 1/a* are neglected 
since a is large. However, under usual conditions of 
heating and cooling, the second term on the right-hand 
side of Eq. (58) can be dropped compared to the first 
term, and therefore / can be further approximated by 


; l 5 i 2A? — B | 
j= — ex sin : 
a(1+B—A)*°P[4°™ (B+ 449)” 
(59) 
Substitution of Eq. (59) into Eq. (56) yields 
Dyytt « 
»= Fs Fe! 
te K(1 +B — A?)'” 
d { 14 (s , 2A? — B 1 
Jexp sin~ 
dx "PLA (B? + 442)? 
sin! a )}t (60) 
(B2 + 4A?)'?/ Hf 
or, upon rearrangement, 
we « D l 
as ae = —-* =e 
bw K* (1+ B — A?’)” 
a’d ™ : (sin~! a + sin! ~) (61) 
tu” f 
where 
2A? — B 1B B 
«= (Bt + 4A” (B? + 442)? 


It is now assumed that the x-variation of the wall tem- 
perature is small compared to the x-variation of the 
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wall friction. 
due to radiation because the heat emission due to 


This is justified in the case of cooling 


radiation varies as the fourth power of the wall tem. 
perature, whereas, in comparison, the heat transfer 
from the fluid [see Eq. (40)] varies directly as the wall 
shear. Certainly, when the boundary layer is insy. 
lated, the wall temperature is constant along the plate, 
If the shear at the leading edge (x = 0) of the plate is 
assumed infinite, the lower limit on a will be zero. 
Hence, integrating Eq. (61) over x and a and assuming 
that the upper limit of a is large, there results 


Pull oX 
Kw 
a D K se | 
i’ ae (sin c sin (62) 
int aR ius6hC(UPC 
or 
oa .. 
7, (sin—! a + sin! B) = const. + 
A (Cyy) 
| ma 
x Crm Rw Tr. (63) 
where 


Eq. (63) can be transformed into an equation in 
terms of the free-stream rather than wall conditions. 
In the first place, since py = p.(7T./T), it follows that 
Ciw = €t., (T/T). Secondly, a power law for viscos- 
ity can be assumed (namely, i = molt. ry 
whence R, = R.(T,,/T.)~"' +”. Eq. (63) then be- 
comes 


ij, (sin~' a + sin~' 8) = const. + 


Ate." (T/T .) 


l l 2 rw 
K (in i> 2 ° in =) (64) 


The constant in Eq. (64) can be determined from the 
requirement that, when M.. = 0 and 7,,/T.. = 1, Eq. 
(64) must reduce to the incompressible fluid case’— 
namely, 

(65) 


1/(c,)'? = 1.70 + 4.15 logw Racy, 


for local skin friction. Adjustment of the constant 
yields the following final equation for local skin-friction 
coefficient when wall temperature, free-stream Reynolds 
Number, and Mach Number are arbitrary: 


0.242 : . 
aan —_——— (sin-! a + sin=! B) = 0.41 + 
A(¢y,.)’(Tw/T a)” 
1 + 2 ‘oo ne 
logio Ra Cr, — - log 10 1 (66) 
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EQUATIONS 66 ano 45 
Tw = 1.0 
To 


Twig 
Tw 


TURBULENT HEAT-TRANSFER COEFFICIENT, Chew 


S 
- * 
o 


REYNOLDS NUMBER, Rao 
Turbulent heat-transfer coefficient for air at Mao = 2. 


Fic. 2. 


The local heat transfer can now be obtained by using 
Eq. (45). 

Figs. 1 to 7 present values of the turbulent heat- 
transfer coefficient computed from Eqs. (66) and (45) 
for Mach Numbers 0, 2, 4, 5, 6, 8, and 10 and various 
In these figures, the highest value 
Values 


temperature ratios. 
of T/T... corresponds to the insulated case. 
used for w and y are 0.76 and 1.400, respectively. 

A formula for the total (mean) skin-friction coefficient 
can also be derived. In terms of wall condition, the 
mean friction coefficient is defined by 


l x 
> CPt ot = f Tx 
2 0 


Direct integration of Eq. (60) then yields, for large 
values of a, 


1 r.\° D a 
9 CR _ (7) K exp F 


Comparison of Eqs. (67) and (62) requires that C;,, > 
Cy, at least for large values of a—that is, for small values 
Of Cm. Bearing this in mind and knowing that the local 
and mean frictions are different at moderately large 
friction coefficients, the form of the equation for 


(sin-! a@ + sin | 
(67) 
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Cy is maintained the same as that of Eq. (63), although 
the constant is different. Hence, 


9'/2 


= -onst. 
A(G,) const. + 


Te (sin—'a@ + sin-'8) = 


aL CroR (7) | (68) 
K N Cp tN T. ) 


or, in terms of free-stream conditions and with the use 
of the power viscosity law, 
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(sin~! a + sin-!8) = const. + 


] 2 Tw» 
K (in R.C;,, — = “in 7) (69) 


The constant in Eq. (69) can be determined, in a man- 
ner similar to that above, when the requirement is 
made that Eq. (69) must reduce to the von Karman 
incompressible, mean skin-friction law—namely, 


0.242/C;,,’* = logiwR.C;., (70) 


when VM, = O and 7,/7T. = 1. Hence, the final 


formula for mean skin friction is 


0.242 (sin—! a + sin-18) 
air Var Pr Ve an : - 
A(Cy,) (Tw/ To)’ 


-T 


* logs a (71) 


‘+ 2 
‘ |e 


logioR co Cr. — 


Figs. 8 to 14 present values of the turbulent mean skin- 
friction coefficient obtained from Eq. (71) for Mach 
Numbers 0, 2, 4, 5, 6, 8, and 10 and various temperature 
ratios. In these figures, the highest values of 7,,/T. 
correspond to the insulated case. The values used for 
w and y are 0.76 and 1.400, respectively. 
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As a matter of interest, the temperature through the 
boundary layer has been computed for M., = 5 and 
R.. = 10° for various wall temperatures using Eqs. 
(34) and (54) in the form 
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TURBULENT BOUNDARY 
where 6 is the boundary-layer thickness and K is taken 
gs 0.40. When the laminar sublayer is disregarded, 
the temperature distribution appears as shown in 
Figs. 15 and 16. 

A relationship between the local and mean friction 
coeflicients can be found in the following manner: 


1 ia | = 
CipPwlto?*X = J Twix = 5 Proll « 2 cr dx 
2 0 2 0 


it follows that 
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1 (Cy 1C fy 
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dx dx (73) 


Differentiation of Eq. (68) with respect to x gives 
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or, upon multiplication of both sides by x and expansion 
of the fight-hand side, 
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But Ry = pullaX/My, Whence x(dR,,/dx) R,,, so that 


Eq. (74) yields, after collection of terms, 
dh 
7——= 
dx 
V2 K , rl 
— | pg (sin7t @ + sin! B) Cy * + Cry | 


: A 
This expression substituted into Eq. (73) gives 
0.558 (1/A) (sin~! a + sin=! B)C;,, 
g= — ; — : = 
; 0.558 (1/A) (sin—! a + sin—18) + 2C,,, ”’ 


where 1/2K = 0.558. Again, C,, = C,,, (T/T) and 
Ow = C7,(T»/T.) so that Eq. (75) becomes 


C—O ——— a - 
 .——1 AH — Sees EQUATION 7| Sint 
si aS SE EY | Pry lero “OT 
s } Yt Te —+++4 

s 14 1} 

r] 

a 

1 | 








HH 

/ } 

1] | 

+ 4 
|| | 
1] | 

[HIE L | 


5, 





TOTAL TURBULENT SKIN-FRIC TION COEFFICIENT, 


10 10” 10° 10 
REYNOLDS NUMBER, Ro 


Fic. 11. Total skin-friction coefficient for air at M, = 5. 
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0.558 (1/A) (sin-! a + sin’ 8)C;,, 


Cha 7 \/s 
eee 1 ° 1 ° 1 ~ Ws i» 
0.558 (sin—! a + sin! B) + 2C,, — 
A r% 
(76) 
TURBULENT SKIN FRICTION ON AN INSULATED PLATE 


When the heat transfer to or from the boundary layer 
is zero, the wall temperature is given by 
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Fic. 15. Temperature distribution through the turbulent bound- 
ary layer for air at Ro = 10’ and Mo = 5. 


(Tig ee = ] + [(y _ 1) 2|M..” 


The mean skin-friction formula [Eq. (71)] becomes, 
then, 


0.242 tes oR . 
nv = + = logi RoC, + 
Cy, 
1 + 2w aes 
logio (1 — A”) (77) 
where 


1—d2 = {1 + [(y — 1)/2]M.2}7 


Eq. (77) is plotted in Fig. 18 for w = 0.76 and y = 
1.4. Also plotted are data obtained at. Daingerfield.” 
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It appears that the data are approaching the writer's 
theoretical curves. 

The Sutherland formula for viscosity 


“ -(Z)"Es5 
a fe T4+S 


In this formula 7 is absolute tem. 


namely, 
(78) 


can also be used. 
perature and S is a constant. According to the data 
of Tribus and Boelter® (see Fig. 17), it is found that 
S = 110°K. Use of this formula in Eq. (68) ulti. 
mately yields, for the insulated-plate case, 


0.242 sy, cin" X 
oe - d*)”? = logic Ruals. a 
Cy. : 
? Od? 
logio (1 — A?) (1 ap .) (79) 
where 0 = S/T... 


The ratio of the compressible to incompressible fluid 
skin-friction coefficient, computed from Eq. (77) at 
Reynolds Number of 10’, is plotted in Fig. 19. For the 
purpose of comparing results obtained using the power 
viscosity law and the more accurate Sutherland law, 
this ratio was also computed from Eq. (79) at the same 
Reynolds Number and is plotted in Fig. 19. The free- 
stream temperature was taken as —67.6°F. It is 
seen that the difference is not great, indeed practically 
negligible. Included in the figure is the friction vari- 
ation for laminar boundary layers. 

von Karman’s results for turbulent skin friction are 
also plotted in Fig. 19. These results are obtained by 
substituting wall conditions directly into Eq. (70) to 
produce 


0.242 


“yy 3 (l — 2)’ ‘= logwR. C,. ob w log 10 (l -— d”) 
C; 5 


x 


(S80) 


using the power viscosity law. The reason for the dif- 
ference between the writer's results and those of von 
Karman is simply that von Karman has assumed that 
the wall conditions apply all the way across the bound- 
ary layer, whereas the writer has allowed the density 
to vary with velocity in the case of constancy of energy 
per unit mass throughout the layer. It seems that von 
K4rman’s formula would be optimistic from the point 
of view of drag. 

It is now advisable to reiterate how von Karman’s 
For lack of more gen- 
Busemann's” 


wall procedure was developed. 
eral information, von Karman* used 
discovery (that for a high Mach Number the velocity 
profile in the case of laminar flow is approximately 
linear) and arrived at an approximation for the vari- 
ation of laminar friction with Mach Number (see Fig. 
20). von Karman then said that this result could be 
approximated by merely replacing the free-stream con- 
ditions of density and viscosity by the wall condition 
in the Blasius formula for laminar incompressible flow. 
Hence, given 
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Fic. 18. Mean turbulent skin friction as a function of Reynolds Number and Mach Number for zero heat transfer. 
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which, for w = 0.76, does cross the linear approxima- 
tion at M.. = 3.5 (see Fig. 20). The fact that the 
linear approximation could be reapproximated by Vah 
using wall conditions in the incompressible flow formula ° 2 4 6 8 10 Rey 
gave von Karman the idea that perhaps a first approxi- FREE-STREAM MACH NUMBER, M,, Eq. 
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-—{l— d 
0 Po Ua U 6 
since 


15 ut ii y 
= pent f tS (1 _ \a (?) 
° pe Re a 6 


and 


Eq. (87) is plotted in Fig. 22 at R. = 10’, assuming 
that Eq. (72) represents the velocity profile and using 
the planimeter to obtain the integral. The von Kar- 
man-Tsien laminar thickness variation is likewise 
shown. 

A turbulent compressible fluid boundary-layer thick- 
ness can also be computed using wall conditions in the 
formula for the incompressible fluid thickness 6;. If 
the logarithmic velocity profile is assumed to hold, 6; 
can be obtained in closed form because the integral 
[Eq. (87)] can be readily evaluated when /., = 0; the 
result is 


F 1 0.558 
é& = — Cy-x- , ores (88) 
2 Cr wg ( = (2c; 4 0.558) | 
where, from Eq. (76), 
cr = 0.558C;/(0.558 + 2C,”*) (89) 


Values of C;,, obtained from Eq. (70) when the wall 
Reynolds Numbers are used, are then substituted into 
Eq. (88) and (89) to yield values of 6... The relative 
thickness obtained in this manner is also plotted in 
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Fig. 22. The constant C [defined by Eq. (83) ], which is 


required to make the incompressible-flow formula [Eq. - 
(88)] yield the compressible-flow thickness [given by 
Eq. (87) ],is0.24at M., = 3. 


DISCUSSION OF THE THEORY 
Recovery Factor 


It was assumed in the development of the foregoing 
theory that the turbulence Prandtl Number c)e/x, as 
well as the laminar Prandtl Number c,u/k, was equal to 
unity. As a result of this assumption, it was con- 
cluded that in the case of the insulated plate the total 
energy per unit mass, defined by Cpl + (u?/2), is con- 
stant across the turbulent boundary layer. This is the 
same conclusion reached for fully laminar layers when 
the Prandtl Number is assumed equal to unity. 

Now, in the case of the insulated laminar layer, the 
fact that the Prandtl Number is not unity, but. rather 
about 0.75, results in a total energy distribution simi- 
lar to the one shown in Fig. 23 in which the energy has 
migrated from regions near the wall to regions near the 
free stream. The curve of Fig. 23 has been computed" 
using the Crocco method to solve the differential 
equations. A manifestation of this migration is the 
usual experimentally observed wall temperature of in- 
sulated plates which is lower than the total free-stream 
temperature. This experimental fact is specified by 
the so-called recovery factor, defined by 7 = (Tx,,. — 
Fits a y= T..), in which Toe is equal to 7, {1 + 
[(y — 1)/2])M."}. With laminar boundary layers, the 
recovery factor, is approximately equal to 
(Np) 2A 

However, in the case of the insulated fuily turbulent 
boundary layer, one would expect that the turbulence 
Prandtl Number would also be somewhat different from 
Indeed, a recovery factor different from unity 


Nam.» 


unity. 
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turbulent boundary layer at Wo = 2.8. 


is observed with turbulent layers. It is found!* that 
Nurb, iS approximately given by (Np,,,.)’° Further- 
more, total temperature measurements'* have been 
carried out at the Aerophysics Laboratory of North 
American Aviation, Inc. (see Fig. 24). These meas- 
urements were made in a '/4-in. turbulent boundary 
layer on the wall of the N.A.A. 50 sq.cm. supersonic 
wind tunnel. Fig. 24 definitely shows a pattern similar 
to that of Fig. 23. It appears that the assumption 
of Prandtl Number equal to unity is a better ap- 
proximation for the turbulent case than for the lami- 
nar. 

The error involved in skin-friction calculation be- 
cause of the assumption of Prandtl Number equal to 
unity can be estimated in the following manner: It is 
first noted that in the case of insulated laminar layers 
the temperature rise for Np, = 0.75 is proportional to 
the temperature rise for Np, = 1 at points throughout 
the boundary layer. It is then assumed that in the case 
of insulated turbulent layers the effect of recovery fac- 
tor on the temperature distribution is roughly the same. 
This leads to the conclusion that an effective Mach 
Number rather than the true Mach Number should be 
used in the above-developed friction formula. Since 
Pin T. {1 + nl(y —1)/2]M.33 and mur. ~ 0.9, it 
follows that the effective Mach Number should be about 
0.95 M., all other properties remaining unchanged. 
The effect of recovery factor on C; would be negligible 
as seen in Fig. 19. This is compatible with the same 
conclusion reached for laminar layers (Fig. 19). Be- 
cause of the recovery factor for turbulent layers, the 
equation for heat transfer would become 
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Qo = CuolpPothal « E ++ —— (0.95 M..)? - i 


(90) 


in which, as previously stated, a more exact value of 
Cu than that given by Eq. (45) requires a more exact 
knowledge of the turbulence mechanism. 


Velocity Distribution 


The turbulent velocity distribution computed by the 
foregoing theory [Eq. (72)] is expected to be too full, 
This is concluded from the fact that, in the case of in. 
compressible flow measurements!‘ on flat plates with 
zero pressure gradients, the velocity distribution seems 
to follow a power law rather than a semilogarithmic 
law. The difference between the law for plates with no 
pressure gradient and the semilogarithmic law for pipes 
with pressure gradients is apparently due to the exist- 
ence or nonexistence of the pressure gradient. (Cer- 
tainly, in the case of laminar flow, the parabolic velocity 
distribution in a pipe is fuller than in the Blasius solu- 
tion for the plate.) However, as far as the derivation 
of a law for incompressible fluid skin friction is con- 
cerned, use of the logarithmic law for velocity yields a 
logarithmic friction law that contains the pertinent 
variables in the proper functional form but whose con- 
stants must be determined by experiment. In the 
development of the above theory for compressible 
fluids, the same procedure is assumed valid. 
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First- and Second-Order Theory of 


Supersonic Flow Past Bodies of Revolution’ 


MILTON D. VAN DYKEt 
The RAND Corporation and California Institute of Technology 


SUMMARY 


Methods are studied for improving the existing perturbation 
theories of supersonic flow past bodies of revolution. Applica- 
bility of the theory at high Mach Numbers is emphasized. 

For axial flow, a second-order solution is found which represents 
a considerable improvement over the first-order result. For 
inclined flow, a second-order solution is not feasible except for a 
cone. Comparison with the exact solutions for cones shows that 
the slender-body series expansion causes large inaccuracies in 
both axial and inclined flow. 

The conclusion that first-order theory predicts the flow no 
better than slender-body theory is shown to be erroneous. When 
first-order theory is properly used, making no unnecessary approx 
agreement is found with exact 
The order estimates used to 


imations, greatly improved 
solutions and with experiment. 
justify the approximations are shown to be invalid in most prac- 
tical cases. A “hybrid’’ theory, combining first-order cross flow 
and second-order axial flow, gives further improvement. A 
physical explanation is advanced for the marked superiority of 
first-order theory over the true “‘linearized”’ theory. 

Nonlinearity in lift is shown to result primarily from viscous 
separation of the cross flow along the after portions of a long body. 
The magnitude of the resulting normal force can be estimated 
with reasonable accuracy using two-dimensional viscous sweep- 
back theory. 


SYMBOLS 


= radius at corner 


a = 

b = reference radius (usually base radius) 

A, 8, C = constants determined by boundary conditions 

c = local speed of sound 

Co, = drag coefficient for circular cylinder 

Cu = pitching moment coefficient 

Cy = normal force coefficient 

Cp = pressure coefficient, (pb — po)/(1/2)p.U? 

Cp,, Cp,, Cp, = coefficients in series expansion of Cp 

Cx = axial force coefficient 

E = complete elliptic integral of first kind with 
modulus k? = (1 — ¢)/(1 + #) 

fig = arbitrary functions 

K = complete elliptic integral of second kind with 
modulus k? = (1 — #)/(1+ 24) 

l = length of body 

M = free-stream Mach Number 

Mp = pitching moment about nose 

N = normal force; also (y + 1)M?/2,? 


Presented at the Aerodynamics Session, Eighteenth Annual 
Meeting, I.A.S., New York, January 23-26, 1950. Revised and 
received September, 1950. 

* The author is indebted to Lester L. Cronvich and K. D. 
Bird, of the Applied Physics Laboratory, The Johns Hopkins 
University, who kindly made available, in advance of publica- 
tion, the experimental results of reference23. Grateful acknowl- 
edgment is due Miss Iris Munson, who calculated many of the 
examples presented in this paper. 

t Consulting Aerodynamicist, RAND, and Research Fellow, 
Cal-Tech. 
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= static pressure 


= radial coordinate 


p 

q = local speed of flow 

r 

R radius of meridian curve of body 


(x) = 
t = conical variable, Br/x 


axial, radial, and azimuthal velocity compo- 


~ 
& 


nents 
= free-stream speed 
axial coordinate 


= axial force 

angle of attack 

VM? -1 

adiabatic exponent 

maximum slope of meridian curve 

azimuthal coordinate 

density of fluid 

= first-order perturbation potential 

= first-order perturbation potential for axial flow 
first-order perturbation potential for cross flow 
second-order perturbation potential 

= second-order perturbation potential for axial 


ll 


SSSESTS* LDR RRG 
| ll | uu ud 


flow 
1 = second-order perturbation potential for cross 
flow 
higher terms in series expansion of @ 


exact perturbation potential 


$2, $3, 1, PS = 


Go = 

x0 = second-order correction potential for axial flow 

vo = particular solution of second-order equation 
for axial flow 

yn’ = partial particular solution of second-order 
equation for cross flow 

Q = complete velocity potential 

( Ja = free-stream conditions 

( ja tbe = values ahead of and behind corner 


INTRODUCTION 


 pooge YNIC FLOW PAST BODIES of revolution at pres- 
ent can be calculated either by the numerical 
method of characteristics or by perturbation theory. 
The characteristics method will, in principle, yield any 
desired accuracy at an arbitrarily high Mach Number. 
The method is laborious, however, and its application to 
inclined flow is not yet well developed. 

Perturbation theory is not applicable above the 
speed at which the Mach cone at any point lies parallel 
to the surface of the body. Even below this limit, 
existing perturbation theory is unsatisfactory in several 
respects. For the axial flow, first-order theory gives 
only fair accuracy in practical cases. Many authors 
have adopted the additional simplification of slender- 
body theory, which further reduces the accuracy. Re- 
cently, the slender-body approximation has been ex- 
tended to second-order, but improvement results only 
for extremely slender shapes. 
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Inclined body of revolution. 


More serious objections arise in the case of inclined 
flow. The lift predicted by the usual first-order theory 
is in poor agreement with the exact solution for cones. 
The slender-body result lies much closer in this case, 
but consideration of other shapes indicates that this is 
a coincidence. Again, the slender-body expansion has 
been extended to second-order, but the series diverges 
so rapidly that it is useful only at moderate supersonic 
Mach Numbers. 


The present paper is concerned with improving upon 
the existing perturbation theory, particularly at high 
Mach Numbers. Second-order solutions are investi- 
gated without recourse to the slender-body approxima- 
tion. The most advantageous use of first-order theory 
is also studied. Finally, one of the principal effects of 
viscosity is accounted for. 


THE ITERATION PROCEDURE 


The first-order (or linearized) theory of supersonic 
flow may be regarded as the first step in an iteration 
procedure. Recently, the second step has been in- 
vestigated for certain three-dimensional bodies." ? 
The iteration equations will now be developed for super- 
sonic flow past inclined bodies of revolution. 


The Perturbation Equation 


Consider uniform supersonic flow past a body of 
revolution inclined at angle of attack a. The boundary 
condition at the body is simplified by choosing a cylin- 
drical coordinate system aligned with the body axis. 
Then it is convenient to resolve the free-stream velocity 


SCIENCES—MARCH, 1951 


into axial and cross-flow components, as indicated jp 
Fig. 1. The shape of the body is defined by its (con. 
tinuous) meridian curve r = R(x). 

The body is supposed to be slender, in the sense that 
the slope of the meridian curve R’(x) is nowhere greater 
than a number «¢ that is small compared with unity, 
(Later, of course, one tries to apply the results to as 
large € as possible.) Then the body is pointed or has 
sharp lip; in the latter case, the internal flow is sup. 
posed to enter at supersonic speed. Under these condi. 
tions, at a given Mach Number, the external flow js 
entirely supersonic for sufficiently small e. 

The flow is assumed to be isentropic throughout, so 
that there exists a velocity potential Q(x, 7, 0). Light- 
hill’ has indicated that this assumption is valid (at 
least at moderate Mach Numbers) to within the ac- 
curacy sought here. The equation of motion is then 


—" , _ %2) 
(c? — 2,°)8,, + (c* — 2,7)0,, + ( c? — — = 
¥ y 
%?\ 2 Q% 2 % 2 
e+ 4 r—-9,—-—-2-—9,— 
yr? , ¢+ r r 
Z02,2., = 0 (la) 


where the speed of sound c is related to c, its value in 
the free stream, by 


Subscripts are used to indicate differentiation. 

A perturbation potential ®(x, 7, @) is introduced in 
the usual way and, for convenience, is normalized 
through division by the free-stream speed 


Q(x, 7,0) = U[xcosa+rsin acos 6 + P(x, 7r,0)] (2a) 


The axial, radial, and azimuthal velocity components 
are then given by 


u/U = 2,/U = cosa+ ®, ) 
v/U =2,/U = sin acos@ + 9, + (2b) 
w/U = Q,/rU = — sin asin 0 + (®,/r) 


Substituting into Eqs. (1) gives the exact perturbation equation. 
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where 6B? = M? — 1. 
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SUPERSONIC FLOW PAST 


Boundary Conditions 


Physical considerations suggest that all perturbations 
should vanish upstream of the body so that 


£(0, 7, 0) = #,(0, r, 0) = 0 (3b) 


and that the flow should be tangent to the surface of 
the body so that 


6,(x, R, 0) + sin acos@ = R’[cos a + ®,(x, R, #)] (8c) 


Actually, the upstream conditions should be supple- 
mented by the shock condition along the bow wave. 
However, just as in the case of plane flow,‘ this refine- 
ment is unnecessary for the order of accuracy sought 


here. 


Pressure Relation and Force Coefficients 

After the potential field has been determined, the flow 
speed g is given in terms of the velocity components of 
Eq. (2b) by 


g?/U? = (u/U)? + (v/U)? + (w/U)? (4) 


Then, because the flow has been assumed isentropic, 
the pressure coefficient is given by 


. 2 4 y¥-1 ( “yy (y-D ) 
= M11 — — —1l} 
2 a U? f 


(9) 

Previous writers have often made approximations in 
the tangency condition, Eq. (3c), and in the pressure 
relation, Eq. (5). The effect of these approximations 
will be considered later; meanwhile the exact relations 
will be used. 

It is convenient to resolve the aerodynamic forces 
action upon a body of revolution into an axial force X, 
anormal force N, and a pitching moment 1/, about the 
nose. Referred to a reference cross-sectional area (such 
as the base) of radius } and to the length / of the body, 
the corresponding dimensionless coefficients are given 


by 
X l ! 
Cy = = R(x)R’(x) X 
, (1/2)p.. Unb? rb? f sini 
& Cp(x, R, 0) dd dx (6a) 
0 
N | a 
Cy = — , — — ~ a f R(x) 
(1/2)p.. U*rb? rh? Jo 0 
Cp(x, R, 0) cos 0d@dx (6b) 
Mp l 7 
Cy = Mr : - y xR(x) X 
(1/2)p.. Urb? rb*l Jo 
(6c) 


| : Cp(x, R, 6) cos 6 dé dx 
( 


The First-Order Problem and Its Solution 


The perturbation equation, Eq. (3a), is completely 
equivalent to the original nonlinear equation of motion. 


BODIES OF REVOLUTION 163 


Consequently, simplifying assumptions must be made 
in order to solve it. The well-known theory of von 
Karman and Moore® for axially symmetric flow and of 
Tsien® for inclined flow is based upon the assumption 
that the right-hand side of Eq. (3a) can be neglected, 
leaving the wave equation. Its solution will be termed 
the first-order potential and is denoted by the script 
lower-case letter ¢. 


The first-order problem is, therefore, 


Yrr + (y,/r) + (Yoo r*) _ Bor: = 0 (7a) 
¢g(0, r, 8) = ¢,(0, 7, 2) = O (7b) 
g(x, R, #0) + sin acos@ = R’[cos a + ¢,(x, R,6)] (7c) 


Lighthill* 7 first noted that this is not the linearized 
problem (as it would be had the x-axis been aligned 
with the free-stream direction) because terms linear 
in ® appear in the right-hand side of Eq. (3a). The 
consequences of this distinction will be considered 
later. 

The first-order problem is satisfied exactly by set- 
ting 


g(x, r, 0) = g(x, 7) cos a + g(x, r) sin acos@ (8) 


The first term corresponds to the axial component of 
free-stream velocity (Fig. 1); the second, to the cross- 
flow component. Previous writers have here assumed 
that the angle of attack is infinitesimal, so that (cos a, 
sin a) may be replaced by (1, a). The effects of this 
approximation will be considered later. 


The first-order problem is thus separated into two 


completely independent problems. For the axial 
flow, 
orr + (¢0,/7) — B*go,2 = 0 (9a) 
¢o(0, 7) = go,(0, 7) = O (9b) 
go,(x, R) = R’[1 + go,(x, R)] (9c) 
and for the cross flow, 
¢irr + (¢1,/r) — (¢1/r?) — Bur, = 0 (10a) 
gi(0, 7) = gi, (0, 7) = O (10b) 
1 + ¢1,(x, R) = R’gui,(x, R) (10c) 


The general solutions of Eqs. (9a) and (10a) which 
also satisfy the upstream conditions of Eqs. (9b) and 
(10b) are® ° 


*x-Br f(z dt 
g(x, r) = — / f(g) a (lla) 
J0 \ (x — £)? i B?r? 
1 fe (x — g(t) dé 
gi(x, 7) = f - (11b) 
er Jo Vx— — pF 


The corresponding first derivatives are 











f'(&) dé conditions, Eqs. (9c) and (10c). This leads to Volterr, 
Por = s wl le = t)? — Br? | integral equations, which can, in general, only be solved 
| 

( 


numerically. The procedure introduced by von Kér. 
man and Moore® and by Tsien® involves superposition 


1 f 
$0, 
r ) 


*x-Br 

“Br (x — £) f"(é) dé 
‘ : 

x-Br 


(x — §)* — f° ( lle of the solution for a cone. The procedure is clearly 
l (x — £) g’(é) dé —— described in the standard texts of Sauer*® and Fer! 

iz = Br J V(x — §)? — BY? | and will not be repeated here. 
l “x-Br (x — £)? g(t) dé | von Karman" introduced an asymptotic solution of 
_— Br? I Wis = £)? — BY? | the integral equations for extremely slender bodies, 
This slender-body theory has been elaborated by British 


The functions f and g are determined from the tangency _ writers using operational calculus.* * 1! !? 


The Second-Order Problem 

It is natural to attempt to improve the first-order solution by iteration. If the first-order solution in the form of 
Eq. (8) is substituted into the right-hand side of Eq. (3a), a nonhomogeneous wave equation results. Its solution 
will be termed the second-order potential and is denoted by the formal lower-case letter ¢. 

This second-order equation, together with the boundary conditions of Eqs. (3b) and (3c), can be satisfied exactly 
by setting 


o(x, r, 0) = go cos a + ¢; Sin a cos 8 + (d2 + ¢3 cos 20) sin® a cos a + (¢4 cos 6 + ¢5 cos 36) sin’ a (12) 


where, as in Eq. (8), the ¢, are functions only of x andr. Again, the problem separates into independent problems 
for each of the ¢p. 

For the present, consider only the first two terms. Making use of the fact that go satisfies Eq. (9a), the second- 
order equation for the axial flow becomes 


po, P 
Pore + ' pa B° doz, = 


: eo - 
2 u:) Poz2P0:" | 


(13) 


. Pe det ee : 
M? 1 2erert + [2 + (y= 1).M?]¢0,,¢02 + P0rrP0r~ + M*Q0,290;" + 207790790; + (1 + 


Similarly, the second-order equation for the cross flow becomes 


oir - > 
Pirr + r -e - B’ diz, = M° 121 + gor) (1 + $17) + 2¢0,,¢0,(1 = ¢1;) + 2¢0,(1 + $02) Piz + $0, *Pirr + 


: ied ‘ > 
M* Q122(2G0, + gor? + Gor?) + Por (2 + Goze) Prez + (CY — 1)M*Qozz[(1 + Goze) G12 + o-(1 + gi) ] + 


{ 
24077( 1 + $0r) Piz + 20021 P0,Piz f (14) 


Series for Pressure and Forces 
In first- or second-order theory, the expressions for pressure and forces can be simplified. Substituting Eq. (12) 
into Eqs. (2b), (4), and (5) and expanding in series gives 


Cp = Cp, + Cp,a cos 0 + (Cp, + Cp, cos 20)a? +... (15a) 


where, on the surface of the body 


2 i =] 
Cm = 45 (O""-"-1), Q=14 T—— M*[1 — (1 + R14)(1 + bo,)?] 
Tm a 


Cp, = —2(1 + R’*)(1 + op) 61,07 ~ ? 
‘\" \ r 
a =. {(1 7 *.) — (1 + R’*) [2(1 + ¢o,)? — oi.” — 4(1 + 6.) x.) i (15b) 
(1/2) M(1 + R’)2(1 + doy) *gie2O® ~— /- 
Coy = (1/2)4 [1 + (¢:/R)]? — (1 + R’*) [br.? + 4(1 + oe) bae]} QO — ? + | 
| 


(1/2) M*(1 + RY %)7(1 + gag)*big2°- PM? 
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Then, according to Eqs. (6), 


c “x _ 


“1 
Cy = (a/b?) | R(x)[ — Cp,(x)] dx + ... 
0 
l 


Cu = (a/b?) f xR(x)Cp,(x) dx + 
0 


ESTIMATION OF ERRORS 


Before considering the second-order solutions, it is 
necessary to discuss the estimation of errors. This 
matter will later prove important in analyzing the effects 
of approximations in first-order theory. 


Slender Bodies at Moderate Mach Numbers 


Actual bounds on the errors in first- and second-order 
theory have not been found. One must be content with 
mere order estimates. British writers have made ex- 
tensive use of order estimates derived from slender- 
body theory. 

If the first-order solution is known, the corresponding 
slender-body solution is the leading term in the series 


expansion in body thickness ratio. For example, on 


the surface of a cone of semivertex angle tan! e, first- 
order theory gives [Eq. (23) ] 
e’ sech—! £8 « ™ 
== (17a) 


$ 
E 


V1i- Be + e*sech™! Be 


Expanding in powers of € and log ¢ gives the slender-body 
result 


gor = —€* log (2/8 €) + Ole! log? €) (17b) 


= O(e log «) 


For any smooth body go, is O(e? log €),* with corre- 
sponding estimates for the other derivatives. (If the 
body has a corner, the flow is locally plane, and all esti- 
mates must be replaced by their two-dimensional coun- 
terparts. ) 

The above expansion implies that fe, as well as «, 
is small compared with unity. For bodies of any 
appreciable thickness, this can be true only if 8 is of 
order unity (Fig. 2a). Hence, the estimates can be 
relied upon only at moderate Mach Numbers. 


Slender Bodies at Hypersonic Mach Numbers 


The iteration procedure breaks down at Mach Num- 
bers greater than that for which Be = 1. Below this 
limit, however, one expects good accuracy fer slender 
bodies from second- and even first-order theory at large 
Mach Numbers (say, such that Be = 1/2). 

O(1), which is the usual defini- 
Both 6 and M are O(e~). 


In this event, Be = 
tion of hypersonic flow.'* 
* The notation y = O(e) is used in the conventional sense that 
lim (y/e) < 
‘> 0 


PAST 


l 
(2/b?) f R(x)R’(x)[Cp.(x) + Cp,(x)a?] dx + 
0 
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(16a) 
(16b) 


(16c) 


The body lies only slightly inside the Mach cone, as 
indicated in Fig. 2b. 
The previous order estimates are no longer applic- 





























able. Thus, on the surface of a cone, first-order theory 
gives 
e? sech—! Be 
SY, = - 7 + O(e*) 
V1 — Be? 
= O(e*) (17c) 
7 
4 
4 
4 
—— -——>-- 
‘\ 
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E<< | 
= O(1) 
Bé<< | } pret 
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=O 
Be = ot) \ B-0(¢) 
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B=0(1) 
Bée =O) 
Fic. 2. Flow régimes for body of revolution. (a, top). Slender 


body, moderate J/. (b, center). Slender body, hypersonic M 


(c, bottom). Thick body. 
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rather than O(e? log €). It is important to notice that 
V1 — p%2 and sech~! Be can no longer be replaced by 
one or more terms of their series developments, be- 
cause the remainder would not be of smaller order. 


Moderately Thick Bodies 


For slender bodies, two entirely separate régimes 
have been distinguished in which different order esti- 
mates apply. Formally, these estimates indicate only 
the limiting nature of the flow as the body thickness 
tends to zero (with V/ held constant in the first case and 
Be held constant in the second). Nevertheless, one 
expects that they furnish physically useful results for 
extremely slender (but nonvanishing) bodies. 

For bodies of moderate thickness, however, the esti- 
mates must be employed with caution. If perturbation 
theory is applicable to a thick body (Fig. 2c), the Mach 
Number is necessarily moderate. At the same time, 
because Be = 0(1), the flow may be termed hypersonic. 
Hence, the two régimes merge. It seems likely that 
order estimates will then prove useful only if one recog- 
nizes that the problem has simultaneously the nature 
of both regimes. Because the hypersonic estimates 
are the more severe, they would govern the problem of 
a thick body. 


SECOND-ORDER SOLUTION FOR AXIAL FLOW 


The Simplified Iteration Equation 


The second-order equation for the axial flow can now 
be simplified. For smooth bodies at moderate speeds, 
the first three terms in the right-hand side of Eq. (13) 
contribute to ¢, quantities of order e' log? «, e* log e, 
and e, while the remaining three terms contribute at 
most O(e® log’ «). At hypersonic speeds, the first four 
terms contribute O(e?); the last two, O(e*). In either 
case the contribution of the last two terms is of the 
order of that that would result from a further iteration 
and can therefore be neglected. The same conclusions 
hold if the body has corners. The iteration equation 
then reduces to 


The velocity components are given by 


Yor = M?*[¢027(¢0 + Nr¢o,) + ¢0;(¢0; + Nr¢oz,) 
M?} go,,(go + Nr¢o,) “ go, [(V + 1) go, + Nr gor, | eis 


Because go vanishes ahead of the body, y will also, and therefore so must xo. 


Yo, = 
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o;; a (do, r)— B° Dorr = 
M*\ 20,02; + [2 = a. Ae 1) M? |g, 0; + 


P07 POrr + l(y = 1)/2] M*g0,7g0,;} (18a) 

with boundary conditions 
o(0, r) = goz(0, r) sid 0 (18b) 
o,(X, R) - R’ {1 Si go, (x, R)] (18¢) 


It is sometimes convenient to regard @p not as the fyl] 
second-order solution but as a correction to be added to 
the first-order solution gp. In this case the iteration 
equation is unchanged, but the tangency condition jg 
modified in an obvious fashion. 


The Particular Solution 


It is necessary only to find some particular solution 
of Eq. (18a). The boundary conditions can then be 
corrected by adding a solution of the homogeneous 
equation, Eq. (9a). Thus, finding any particular solu. 
tion reduces the second-order problem to an equivalent 
first-order problem. 


The particular solution for the axial flow will be de- 
noted by Yo and the correction potential by xo, so that 


go = Wo + Xo (19) 

It was shown in reference 1 that an approximate 

particular solution of Eq. (18a) is given in terms of the 
first-order sclution by 


Yo = M?[¢0,(¢0 + Nrgo,) aad (1/4)r¢o,*] (20a) 


where 


N = [(y + 1)/2](M?/6") (206) 


At moderate speeds, this function satisfies the itera- 
tion equation except for terms of the order of those 
At hypersonic speeds, it fails to 
However, 


already neglected. 
account completely for the triple products. 
examples show that the actual magnitude of the error 
in this case is small and that a complete particular 
solution cannot be expressed in terms of the first-order 
solution, so that the approximate solution is acceptable. 


— (3 4)rgo,,¢0," | (20c) 


(1 4) g0,?(¢o, + 3r¢o,,) | (20d) 


The problem for xo is then just the 


first-order problem, Eqs. (9), with the tangency condition replaced by 


Yo, (x, R) a Xo(X, R) = R’[l + Yor(X, R) oe Xo0,(X, R)] 


Solution for a Cone 


A cone is one of the few simple shapes for which the first-order solution can be written in closed form. 


venient to introduce the conical variable 


(21) 


It is con- 
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t = Br/x (22) 
Then, ¢ = 1 on the Mach cone, and on the surface of the cone ¢ = Be, where ¢ is the tangent of the 


in place of r. 
The first-order potential is known to be® ® 


semivertex angle (Fig. 3). 


yo/x = —A (sech t — V1 — 2) (23a) 


where 
2 


r : (23b) 
V1 — Bre? + €? sech—! Be 


A= 


The particular solution Yo is then given by Eq. (20a). To this must be added the correction potential xo, which is 
some multiple of the first-order solution. Thus, the second-order perturbation potential and its derivatives are 


found to be! 


/ / 1,,a-#)” 
wa —C(sech-'t — V1 — #?) + A?M? | (sech 14)? — (N+ 1) V1 —#sech"?t — gh A 2 | (24a) 
. 2 


ce _ . ; sech~! f 1 3... V1l-# 
go, = —Csech—'t + A?M? | (sech— #)? — (N — 1) — —(N+1)- 4 B°A - (24b) 
V1i-# , 
Wi-# V1 — fsech-'t t sech—" t 
do, = C 4 aur| - 2 = 4iWwen-4+@~-9— 
B t t t V1—-# 


l 9 9 2 V1 aa e 9Ap 
; PAZ +) —, (24c) 


The constant C is best determined numer'’cally from the tangency condition, Eq. (18c), rather than from the 


cumbersome expression that could be written for it. 





30 : ——_—_—_— 


| | 
------ FIRST ORDER | | 
| ——— SECOND ORDER | | 


EXACT 


Second-order values of surface pressure coefficient 
are compared in Fig. 4 with the first-order values and | 
with the exact values!‘ for cones of 10° and 20° semi- 
vertex angles. The second-order solution is seen to | 





represent considerable improvement over the first- 
order solution throughout the range of Mach Numbers 


up to the hypersonic limit. 





Method of Solution for Smooth Bodies 





In general, the first-order solution can be determined 
only numerically. The method introduced by von 
Karman and Moore® involves step-by-step superposi- 
tion of the solution for a cone. To second-order, 
however, the cone solution would produce false pressure 
jumps on a smooth body so that a smoother basic 








solution must be used. 
The procedure will be outlined first for a smooth 
pointed body. The axis of the body is divided into 
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‘\ Fic. 4. First- and second-order solutions for pressure on a cone. 
(a, top). 10° cone. (b, bottom). 20° cone 














Fic. 5. Method of solution for smooth pointed body. 











Body with a corner. 


Fic. 6. 


intervals by choosing points with abscissas Z, as 
shown in Fig. 5. From each point, the Mach line is 
drawn intersecting the body at the point P,. 

The first-order solution is started with the potential 


for the tangent cone with origin Zp) at the nose. This 
potential and the required derivatives are 

go = —Ax(sech-?t — V1 — #) | 

gor = — Asech't | 

yo, = BA(V1 — £/t) ; 

a: (25) 

Yor = — (A/x)(1/V1 — 2°) 

G0rr = (BA/x)(1/t V1i— F) | 

yo,, = — (8A/x)(1/2V1 — #) | 


where A is given by Eq. (23b). All these quantities 
are calculated at each of the points P,. 

Next, the smooth basic solution is added with ori- 
gin Z, at the nose. This potential and its derivatives 
are 





0 = 
r 
20 ie 2t i 
Yor al K 
7 r L+tf 
7 28C a 2t (* +i 
sad T , 1l+t t 
C1 a 7 2t 
rz ( 
” rx Wri1—¢ Vit+e 
BC 1 a | 2t (; 
P0rr 
w x rl-t 1+itv\t 
BC 1 a l 2t 
P0rr _ 
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Dea . 30/ -] } 
— Cx? 1+ = t? } sech!t — — 1-f 
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Yor, = —2C sech' t 

Yor = 2BC(V1 — ?/t) 

Yor = —B°C((V1 — #2/t2) + sech- ft] 


The constant C is chosen so that at point P» the com. 
bined contributions of both solutions satisfy the tan. 
gency condition, Eq. (9c). 

The solution is continued by adding this basic solution 
with origin shifted to Z2, Z3, etc., and determining the 
strength C in each case from the tangency condition at 
the points P;, P,, etc. Adding the contributions of al] 
solutions gives, at each point P,, the first-order po- 
tential g and its derivatives. 

The velocity increments due to the particular second- 
order solution yp can then be calculated at each point 
from Eqs. (20c) and (20d). 

Finally, the correction potential xo and its first de- 
rivatives are determined by repeating the procedure 
used to find go, with new constants C chosen so as to 
satisfy the second-order tangency condition, Eq. (21), 
at each point P,. 

Adding the contributions of Yo and xo gives the com- 
plete second-order perturbation velocity components 
at each point. The pressure coefficient can be then 
calculated from Eqs. (4) and (5). 

The computing labor required for a second-order 
solution is several times that for a careful first-order 
solution. 


Treatment of Bodies with Corners 


Suppose that the body has a discontinuity in curva- 
ture at P,, as indicated in Fig. 6. The procedure just 
described for calculating the first-order solution gives 
poor accuracy behind the corner unless the intervals 
are chosen extremely small. However, the corner cam 
be effectively removed by adding, with origin at Z,, 
the following corner solution: 


x) ~ ac | 
. r (27a) 
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Here, K and £ are the complete elliptic integrals of the 
first and second kind, with modulus k? = (1 — #)/(1 + 
)), The approximate expressions give the values just 
behind the corner. If the reference length a is taken 
to be the radius of the body at the corner, the tangency 
condition is satisfied by putting 


(R,’ age R,') {1 + (goz)a] 
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Here, R,’ and R,’ are the slopes just ahead of and be- 
hind the corner, and (g,)a is the value ahead of the 
corner. The first-order solution can then be continued 
as before. 

In the second-order solution, Yo will be found to be 
discontinuous at the corner, and the solution is sub- 
sequently incorrect. This discontinuity can be can- 
celed by including in xo the following step solution with 


C= - (27b) a : 

B+R, origin at Z;: 

xy = (2C/n) Va/r V2t/(1 +8) K oe. 
C1 2 2t Cc 

= K-E = — — 

sae wx rl-—t 1+? \ ) 86r | (28) 
C1 l 2 /1 ; 3 

x» = —§ - ((e-K)~-< 
wx rl1-t 1+tv\ 8r 


From the approximate forms valid just behind the 
corner it is seen that C must be equal and opposite to 
the jump in Yo at the corner. Also, xo will include a 
corner solution, Eq. (27), with a new constant. 

The external flow past open-nosed bodies can be 
treated, provided the internal flow enters at supersonic 
speed, by starting with the corner solution. 


Computational Details 


As the solution proceeds far along a body, the desired 
quantities are found as relatively small differences of 
increasingly large numbers. In order to obtain results 
valid to three significant figures, it is therefore neces- 
sary, in typical problems, to work with five or six sig- 
nificant figures. 

Simple rules can be given for choosing the intervals 
along the axis. Dimensional reasoning indicates that, 
in order to achieve a given degree of accuracy, the first 
interval for a smooth ogive should be proportional to 
p/8, where p is the radius of curvature of the meridian 
curve at the tip. Experience has shown that an ac- 
curate solution is obtained at all Mach Numbers if the 
first interval is taken as 0.025 p/8. Similarly, if corners 
are subtracted out, it is sufficient, except near the 
hypersonic limit, to choose subsequent intervals no 
greater than 8/2 times the local radius R. 

Before calculating Yo, it is helpful to check the first- 
order solution go by verifying that Eqs. (9a) and (9c) 
are satisfied numerically at each of the points P,. 

If the body has a corner or a discontinuity in curva- 
ture, plots of do,,, $0;,, and @o,, (and, hence, also of 
Yor, Wor» Xor, and xo,) will oscillate somewhat, but the 
final values of ¢o, and @o, will be smooth. 


Comparison with Characteristics Solutions 


The accuracy of the second-order solution for axial 
flow can be evaluated by comparison with examples 
computed using the method of characteristics. For this 


purpose, the examples summarized in reference 15 are 
useful, as well as those of reference 16. 


The pressure distribution on the 12'/2 cal. A-4 ogive 
was calculated at a Mach Number of 3.24. This repre- 
sents a severe test of the method, since the Mach angle 
is then only 10 per cent greater than the cone angle. 
Intervals were chosen such that the points P, lay at 
0.1, 0.25, 0.5, 1, 2, and 3.5 cals. The result is seen in 
Fig. 7a to agree with several characteristics solutions 
summarized in reference 15 to within the accuracy of 
the characteristics method. The second-order solution 
yields too high a value at the tip (as do the numerical 
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Fic. 7b. Comparison of second-order and characteristic solutions. | | 
10° cone-cylinder at M = 2.0750. | | | | 
re) i J 
1 2 3 “ 
solutions) but appears to overcome this error almost at ‘ 
Fic. 8. Effect of expanding in series upon second-order solution 
once. for cone. (a, top). 10° cone. (b, bottom). 20° cone. 


As a test of the modification required after a corner, 
the second-order solution was also calculated for a 10° Broderick!! obtained a second-order solution for sufi- 
cone-cylinder combination at a Mach Number of 2.075, ciently smooth bodies which is based upon the slender- 
which was treated numerically in reference 16. Again, body approximation. It can be shown that Broderick’s 
the second-order result is seen from Fig. 7b to agree solution also results from expanding the present second- 


with the characteristics solution. order solution in series and keeping two terms. 
a it hina If the solution for the cone, Eqs. (23) and (24), is 
Rpoas af Sapanting te Tevtes expanded in powers of ¢ and log 2/t, if the constant is 


It was noted before that the slender-body solution is evaluated from the tangency condition, and if the pres- 
obtained by expanding the exact first-order solution in sure relation is also expanded in series, the result is 
a series and neglecting all but leading terms. Recently, that on the cone 


2 2\3 2 M4 13 l 
Cp = é (2 log Be _ ) + e! | 36 (tog *) — (5M? — 1) log Be + (y¥ + 1) 32 + 7 M? + | + 


9\3 
O | é | log (29) 
Be 
This agrees with Broderick’s result for the cone. 

This series expansion is compared in Fig. 8 with the original form of the second-order solution. Expanding in 
series is seen to cause great loss of accuracy, so that for moderately wide cones the second-order result is actually 
inferior to the first-order solution (Fig. 4). This tends to confirm the earlier conjecture that thick bodies are 
governed by the hypersonic order estimates; according to this, series expansion is unjustified. 


SECOND-ORDER PROBLEM FOR CRosS FLOW 


The Simplified Iteration Equation 


As in the case of axial flow, the second-order equation for cross flow can be simplified. For smooth bodies, terms 
in the right-hand side of Eq. (14) can be neglected if their contribution to ¢), is no greater than O(e® log ¢) at 
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moderate speeds and O(e*) at hypersonic speeds, because a further iteration would contribute additional terms of 
these orders. In addition, for bodies with corners, all quadratic terms in derivatives of g must be retained. Thus, 





the iteration equation simplifies to 


oi, 
di + ; nse 


1 9 ‘ 2 I 2 Y 
y2 a B*hizy = 2M? (w(t + ¢1,) + $0;P irr a $07, ¢0,(1 . ¢1,) + 9 $0; Pirr + 9 


— |] 
M?* X 


; = 
{Perle + go(1 + gir) + Pure (a. Ts ont) FH Por2Piz + P02Pizr + a (30) 


The last three terms can be omitted for smooth bodies, and the terms multiplied by (y — 1) are required only 


for hypersonic speeds. 


The Partial Particular Solution 


A partial particular solution of Eq. (30) is given in terms of the first-order solution by 


2 


1 M? 
y,’ = M? | re + gogo, + 5, 1(0r? + B*y02”) + (Gopir + Pore) + Nr( gore, + ewe) | (31) 


28 


Thus, the problem is reduced to finding a particular solution of 


oi, 


1 ae. 7 9 l 9 
tot ~~ Ot = 2M? ' P0rr Por Pir F 5 Por Pirr — 


It does not seem possible to express a particular solu- 
tion of this reduced equation, even approximately, in 
terms of the first-order solution. Consequently, a 
second-order solution for the cross flow past an arbitrary 
body is not feasible. 


Solution for a Cone 


The simple example of a cone indicates the nature of 
the second-order solution. For a cone of semivertex 
angle tan—! e (Fig. 3), the first-order cross-flow potential 


is given by® ° 
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77s. a 
N¢o;; = + - M? | exe +- ¢1,) + eee | 
r 2 2 
(32) 


where 


Be® 


- _ : (33b) 
(1 + 2e2) V1 — Be? + Be? sech—! Be 


Because the flow is conical, Eq. (32) is reduced to a 
total differential equation by setting 


oi(x, 7) = xw(t) (34) 
The resulting equation can be solved directly; the 
resulting lengthy expression for ¢, need not be repro- 
duced here. 

In this way, second-order values of normal-force 
slope and surface velocity increment due to cross flow 
have been calculated for a 10° cone and are compared 
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in Fig. 9 with the exact solution of Stone as calculated 
by Kopal’s group." The usual first-order result, as 
given in Ferri’s book,’ is also shown. Although the 
second-order surface velocity is accurate, the normal- 
force slope is too great at high Mach Numbers. This 
discrepancy indicates that the assumption of potential 
flow is inaccurate, as noted by Lotkin." 


Effect of Expanding in Series ‘ 


If the second-order solution for a cone 
in powers of ¢ and log 2/t, it is found that! 


is expanded 


oe) > ‘ (2 p 2 : ) 
d = 2—¢(2M*log— — 3M? +2 
( a € og . + + 
9\2 
O |. (10g =) | (35) 


The approximate value of 2 was given by Tsien;* the 
above result was also obtained by Lighthill,* who 


has developed a second-order theory for inclined 
bodies based upon the slender-body approxima- 
tion. 


The first- and second-order slender-body values for 
normal force and surface velocity are shown in Fig. 9. 
As in the case of axial flow, series expansion of the 
second-order solution is seen to reduce its accuracy 
considerably. 


PROPER USE OF FIRST-ORDER THEORY 


Apparently, a second-order solution for cross flow is 
not feasible. However, it will now be shown that first- 
order theory, if properly used, yields much better results 
for inclined was heretofore considered 
possible. 


bodies than 


Approximations in First-Order Theory 


All previous treatments of first-order theory have 
involved some of the following approximations: (1) 
approximate pressure relation; (2) approximate tan- 
gency condition; and (3) infinitesimal angle of attack. 
Slender-body theory includes all these, together with 
further approximation. The justification advanced for 
making these approximations is that the order of the 
error committed is no greater than that due to simplify- 
ing the equation of motion. Hence, Lighthill'!? and 
Laitone’® conclude that exact solution of the first- 
order equation is fruitless and that first-order theory is 
incapable of yielding results more exact than those of 
slender-body theory. 


These conclusions may be valid for extremely slender 
bodies at moderate Mach Numbers. However, at hy- 
personic speeds or for bodies of reasonable thickness, it 
has been suggested that the order estimates underlying 
these arguments cease to apply. Therefore, it should 
not be surprising to find that the conclusions are in- 
correct in most practical cases. 
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Approximate Pressure Relation 


According to Eq. (5), the exact pressure relation jg 
9 -_ 2 Vy 1) 
mts By | 

~ <i + M?{1 — — <. 
yM? \ 2 U? : 


At moderate Mach Numbers, the term [1 — (¢° U’)] 
is O(€* log €), as is the entire second term in the bracket. 
It is therefore permissible to expand the bracket jn 
series, which yields the approximate pressure relation 


Cp ~ —26, + 4, 


Co = 


(36a 


The corresponding approximate forms of the first two 
terms in the series of Eq. (15) are 


Cry ~ —2do, + $0,” \ 


Cp, = —2¢1, 


Lighthill!* has shown that the terms retained here are 
consistent with the first-order theory. 

At hypersonic Mach Numbers, [1 — (g?/U?)] de. 
creases to O(e”), but J/* increases to O(e~*). Hence, 
The series ex- 


(36b) 


the second term in the bracket is O(1). 
pansions can no longer be justified on the basis of order 
estimates. 


Approximate Tangency Condition 


In the exact tangency condition of Eq. (3c), the term 
®, is small compared with cos a. If it is neglected, the 


approximate tangency condition becomes 


(37a) 


®,(x, R, 0) + sin acos@ = R’ cosa 


For the first-order axial- and cross-flow problems, 
Eqs. (9c) and (10c) become 


$o,(x, R) = R’ \ 
gi,(x, R) = —1] 


(37b) 


These approximations have the advantage that the 
resulting first-order solutions satisfy the supersonic 
similarity law (the supersonic counterpart of the 
Gothert rule”’), which is not true if the exact tangency 
condition is used. Thus, the solution calculated for a 
given body at one Mach Number applies to a whole 
family of similar bodies of different thickness ratios at 
correspondingly different Mach Numbers. 

The error involved in the previous approximation 
was seen to increase with Mach Number so that the 
approximation is unjustified at hypersonic speeds. This 
is not true of the approximate tangency condition. 
Except possibly for extremely thick bodies, the resulting 
error must be small at all Mach Numbers. 


Infinitesimal Angle of Attack 


If the angle of attack is infinitesimal, sin a and cos 
are replaced throughout by a and 1. First-order the- 
ory then predicts only a linear variation of pressure 
with angle of attack. This approximation is unneces- 
sarily restrictive, however, because first-order theory 
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also gives the quadratic variation correctly, at least at 
moderate Mach Numbers—that is, a first approxima- 
tion to Cp, and Cp, in Eqs. (15) is obtained by omitting 


g, and ¢3,- This conclusion is implicit in the work of 


Lighthill.* 


At hypersonic Mach Numbers, omission of ¢2, and 
a, cannot be justified by order estimates. An example 
considered below (Fig. 12) suggests that the omitted 
terms become important (but not predominant) at high 
Mach Numbers. Similarly, the order estimates of 
Bq. (38) indicate that cubic terms in a may assume im- 
portance at hypersonic speeds, but the examples below 
suggest that their effect is small. 


Comparison with Results for Cones and Ogives 


The effects of these three approximations and of the 
sender-body approximation will now be evaluated by 
comparison with inclined cones fer which an exact 
theory and some experimental results are available. 

The theory of Stone” gives the exact initial normal- 
force slope for cones. The first two of the approxima- 
tions listed above affect this value. 

Fig. 10 shows the exact solutions for 10° and 20° 
cones, the exact first-order results, and the first-order 
values involving each of the first two approximations 
listed above. Without approximation, first-order the- 
ory is seen to yield much closer agreement than has 
been found heretofore. The effect of linearizing the 
pressure relation is disastrous, while using the approxi- 
mate tangency condition has slight effect. 

The first- and second-order results of slender-body 
theory are also shown in Fig. 10. The simple value of 
2is accurate for slender cones, but this must be a coin- 
cidence peculiar to cones, because adding a cylinder 
behind the cone contributes considerable additional 
lift that is not predicted by the simple slender-body 
theory. 

The exact value of this additional lift just behind the 
corner of a cone-cylinder combination can be calculated 
by combining Prandtl-Meyer expansion theory with 
Stone’s solution. The exact value for the rate of in- 
crease of normal-force slope behind the corner is com- 
pared with the various first-order results in Fig. 11. 
Because the flow is locally two-dimensional, the agree- 
ment is poorer in all cases than on the cone. Again, 
however, the approximate pressure relation causes loss 
of accuracy, while the approximate tangency condition 
happens to improve the agreement somewhat. 

The first-order slender-body value is zero; the corre- 
sponding second-order result has not been derived. 
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The quadratic terms are important if the angle of 
attack is comparable with the maximum body slope, e, 
as noted by Allen.*!_ This is clear from the order esti- 


mates for the terms in Eq. (15) 


Cp = Cr + Cpacos@ + (Cp, + Cp, cos 20)a? +... ) 
Moderate M: € loge €a a? are (38) 
Hypersonic 1: é €a a? ave! 


The effects of assuming infinitesimal angle of attack 
cannot be evaluated from initial-slope values. Re- 
cently, Stone has extended his exact theory for cones to 
include squared terms in angle of attack, and computa- 
tions have been published by Kopal’s group.”? Fig. 12 
compares first-order values with the exact values for the 
coefficients of a? in the series for surface pressure, Eq. 
(15). First-order theory is seen to be superior to the 
slender-body results,* although the agreement decreases 
near the hypersonic limit. The approximate tangency 
condition has little effect and is not shown. The 
approximate pressure relation yields zero for both 
coefficients, and the same result is implicit in the as- 
sumption of infinitesimal angle of attack. 

It is of interest to know whether including only 
squares of the angle of attack suffices to predict the 
pressures on a highly inclined body. For this reason it 
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is fortunate that pressure measurements over inclined 
cone-cylinder combinations are available from tests 
conducted by Cronvich and Bird.** The measured 
distribution of pressure around a 10° cone at the 
highest angle of attack, 12°, and a Mach Number of 
1.50 is shown in Fig. 13. The first-order solution agrees 
well with the experimental results. Approximating 
to the pressure relation completely destroys the agree- 
ment, while use of the approximate tangency condition 
has little effect. The approximation of replacing (sin 
a, cos a) by (a, 1), which is considered here for the 
first time, is seen to shift the pressure distribution 
without appreciably distorting it. Therefore, this 
approximation impairs only the predicted value of axial 
force, the normal force and pitching moment being un- 
affected. 


The effects of the various approximations are similar 
for shapes other than cones. Fig. 14 shows the experi- 
mentally determined distribution of normal force along 
the 12'/.-cal. A-4 ogive at 10° angle of attack and a 
Mach Number of 3,24. The experimental points and 
faired curve are taken from reference 24. The first- 
order theory lies within the experimental accuracy ex- 
cept near the nose, where it is somewhat too low in ac- 
cordance with Fig. 10. Approximating to the pressure 
relation causes considerable error over the forward por- 
tion of the ogive. On the other hand, the slender-body 
approximation becomes worthless over the rear portion. 
Approximating the tangency condition (which is 
not shown for the sake of clarity) has a negligible 


effect. 
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A Hybrid Theory 


The comparison above was restricted mainly to lifting 
effects; poorer agreement results in the case of axial] 
flow alone (though using the Jinearized pressure relation 
Cp, = —2¢0, happens to compensate for the error gt 
moderate Mach Numbers). It is natural (though 
contrary to the views of previous writers) that first. 
order theory predicts the cross flow more accurately 
than the axial flow because smaller disturbances are 
involved. 

This suggests that it is more important to refine the 
axial flow than the cross flow; fortunately, just this can 
be done using second-order theory. A hybrid theory js 
therefore suggested, involving the first-order cross flow 
combined with the second-order (or other more exact) 
axial flow. 

The results of such a hybrid solution are shown in 
Fig. 15 in comparison with the exact results for cone- 
cylinder combinations. Comparing with Figs. 10 and 
11 shows an improvement over first-order theory. The 
superiority of hybrid theory to second-order theory 
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Fic. 14. Comparison of theoretical and experimental normal- 


force distributions along A-4 ogive at I = 3.24, a = 10°. 


(Fig. 9) at high Mach Numbers must perhaps be re- 
garded as fortuitous. 

As a second test of the method, the pressure distri- 
bution has been calculated for a 10° cone-cylinder 
combination at a Mach Number of 2.00 and angle of 
attack of 12°. The result is compared in Fig. 16 with 
the experimental measurements of reference 23. Again, 
theagreement withexperiment is excellent over the cone. 
Immediately behind the corner, however, on the bot- 
tom, the hybrid solution departs considerably from the 
measured result. Since the changes at a corner are 
locally two-dimensional, the pressure just behind the 
corner has been calculated using a Prandtl-Meyer ex- 
pasion, assuming that the theory is exact on the cone. 
It is seen that this result coincides with the measured 
value. Consequently, a truly nonlinear phenomenon 
arises at this point which the approximate theory cannot 
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account for. However, the discrepancy subsides rapidly 
so that the error in lift would be slight. 

Otherwise, the agreement between theory and ex- 
periment is satisfactory with the exception that along 
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TABLE 1 
Body Coordinates Wind Coordinates 
First-order Cle: = 0 Ole = —2M*¢0,_ 
Linearized 7" = 2M Cle: = 0 


where []¢: = ¢1,, + (¢:,/r) — (¢1 iy?) — Bei, 


the top of the body the measured pressures lie con- 
sistently below the theoretical values. It will be shown 
in the next section that this is apparently the result of 
viscous separation. 


Comparison of First-Order and Linearized Equations 


It was pointed out previously that the first-order 
equation is not the linearized equation, although this 
distinction has been largely ignored. For axial flow the 
two equations are the same, but for cross flow the lin- 
earized equation in body coordinates contains an addi- 
tional term, as shown in Table 1. 

A particular integral that reduces the linearized 
equation to the first-order equation is given by the first 
term of Eq. (31) 


VY = M?rgo, (39) 


Lighthill’ has shown that, for extremely slender 
Mach Numbers, the difference 
and first-order solutions is 
hypersonic speeds, however, 


bodies at moderate 
between the linearized 
of smaller order. At 
and for bodies of practical thickness, his order esti- 
mates are invalid; actually, the difference becomes 
great. 

Invariably, the result of solving the linearized equa- 
tion is far inferior to the first-order result. A single 


example will suffice. Fig. 17 compares the first-order 
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Again, as in Eq. (27a), A and E are complete elliptic 
integrals of modulus k? = (1 — t)/(1 + #4). The approxi- 
mate forms hold just behind the corner. 

It has been seen that the approximate tangency 
condition causes no appreciable loss of accuracy. <Ac- 
cording to this approximation, a corner behind the nose 
has no effect upon the cross flow. Therefore, only in 
the case of an open-nosed body must a unit multiple of 
Eq. (40) be added at the start of the solution. 


and linearized solutions for initial normal-force slope of 
a 10° cone. 

A simple physical explanation can be given for the 
superiority of first-order over linearized theory. As the 
body is rotated through zero angle of attack, the char. 
acteristic surfaces for the linearized equation remain 
aligned with the free stream, while those for the first- 
order equation rotate with the body. In actuality, near 
the body the characteristic surfaces rotate slightly 
faster than the body, although far from the body they 
may, at moderate Mach Numbers, move very little. 
Flow disturbances occur primarily near the body and 
are governed by the characteristic surfaces, so that the 
first-order equation represents the flow more accurately 
than the linearized equation. 

A mathematical explanation is more difficult to give 
on the basis of order estimates alone. It seems clear, 
however, that the extra term in the linearized equation 
is counterbalanced by nonlinear terms in Eq. (30). 
For example, using the approximate tangency condi- 
tion, the terms 21/79 ,,(1 + ¢1,) vanish on the surface 
of the body. 


Corner Solution for Cross Flow 


In the step-by-step calculation of first-order cross 
flow, large steps may be used if, as in the case of axial 


flow, corners are effectively removed. A_ suitable 
solution for this purpose is given by 

= tk) = (0 
: t Kk) _ : (40) 
1+t B | 


2-1 
EK - K|~-1 } 
2(1 + 2?) 


Viscous EFFECTS 
Nonlinearity in angle of attack has so far been largely 
excluded by ignoring terms beyond ¢; in the series for 


¢@, Eq. (12). 
in the slender-body approximation by Lighthill,’ who 


The higher terms have been considered 
finds that their effect is slight. The experimental re- 
sults of Figs. 13 to 15 also suggest that the higher terms 
are negligible. 
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SUPERSONIC FLOW PAST 


Most experiments ‘show a lift-curve slope that in- 
creases With angle of attack. This nonlinearity must 
result from effects. not accounted for in the potential 
equation—in particular, from viscosity. An attempt to 
include viscous terms in the equation of motion would 
lead to formidable difficulties. However, it appears 
that one of the most important effects of viscosity can 
be analyzed approximately with good accuracy. 


Lift Due to Viscous Drag of Cross Flow 


Consider an inclined body of revolution consisting of 
a nose followed by a long cylinder (Fig. 18). Far down- 
stream of the nose, the flow may be expected to ap- 
proach that past an infinite circular cylinder. Jones* 
has shown that, for a yawed infinite cylinder, Buse- 
mann’s sweepback concept applies to the viscous, as 
well as the potential, flow—that is, the axial velocity 
component can be disregarded, the essential features 
of the flow being determined by the cross-flow com- 
ponent of velocity. Hence, on the downstream portion 
of the body, flow separation will occur over the upper 
half, corresponding to that over the rearward half of a 
circular cylinder. The resulting force N (Fig. 18), 
which is the drag of the cross flow, appears as an addi- 
tional normal force on the body. 

The magnitude of this force can be calculated from 
the value of the drag coefficient, Cp, for a circular 
cylinder. Taking account of the differences in refer- 
ence area and reference velocity for Cp, and Cy, it is 
found that the separation contributes an increment in 
normal-force coefficient per unit of body length 


given by 


d 2 R 
(ACy) = - Cy = sin’? a (41) 
dx T b? 

where } is the reference radius for the body. The value 
of Cp, to be used here is that appropriate to the effec- 
tive values of Mach and Reynolds Numbers. These are 
the free-stream values reduced by the factor sin a. 
However, for moderate angles of attack and low super- 
sonic Mach Numbers, the effective Mach Number will 
be less than 0.4 (the critical value for a circular cylinder) 
below which Cp, is little affected by compressibility. 
Moreover, Cp, is nearly constant over a wide range of 
Reynolds Numbers. Hence, in most cases, Cp, can be 
taken to be from 1.0 to 1.2. At high angles of attack 
and large Reynolds Numbers, the effective Reynolds 
Number may exceed the critical value for a cylinder 
(approximately 2 X 10° based on the diameter), and 
these values must be modified. 

The above discussion refers to a section of the body 
far down the cylindrical portion; the question arises as 
to where along the body separation first becomes im- 
portant. An approximate answer appears to be that the 
cross flow separates when it encounters an adverse gra- 
dient comparable with that around a circular cylinder. 
Thus, Fig. 13 indicates that no appreciable adverse 
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Fic. 17. Comparison of first-order and linearized solutions for 
normal force on 10° cone. 
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Fic. 18. Cross flow past a long inclined body. 
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Fic. 19. Comparison of theoretical and experimental values 
for nonlinear part of normal force on cone-cylinder combinations 
at M = 2.00 
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pressure gradient occurs around a cone even at relatively 
high angles of attack. Consequently, the flow remains 
attached, as is indicated by the agreement between 
theory and experiment on the leeward face of the cone 
(Figs. 13 and 15). However, if the cone is followed 
by a cylinder, the gradient becomes strongly adverse 
close beyond the corner, and Fig. 15 suggests that 
separation begins almost at once. 

Thus it seems likely that the theoretical pressure dis- 
tribution calculated according to the methods recom- 
mended previously will provide a rough criterion for 
the onset of viscous separation. 


Cross-flow separation has also been investigated by 
Allen.”! Similar concepts were applied to wings of low 
aspect ratio by Kriesis.*® 


Comparison with Experiment 


For the cone-cylinder combinations discussed above, ”* 
experimental values of the nonlinear part of the normal 
force were determined by subtracting the values corre- 
sponding to the initial slope of the normal-force curve 
(which were close to those predicted by either first- 
order theory or its hybrid refinement). These are 
compared in Fig. 19 with the values calculated from Eq. 
(41), assuming that separation starts at the corner with 
full strength. The satisfactory agreement suggests 
that nonlinearity in normal force results primarily from 
viscous separation and that its effects can be predicted 
with reasonable accuracy on the basis of two-dimen- 
sional viscous sweepback theory. 


CONCLUDING REMARKS 


It has been shown that, contrary to the common 
view, supersonic perturbation theory predicts the flow 
past inclined bodies of revolution with good accuracy. 
First-order theory, if properly used, is applicable at high 
Mach Numbers and for relatively thick bodies. The 
wave equation in body axes must be solved, rather than 
the true linearized equation. The usual approximations 
to the pressure relation cannot be made. If pressure 
distributions are required, the angle of attack should 
not be assumed infinitesimal. 

For bodies of other than circular shape, it seems likely 
that comparable accuracy would result from solving 
not the linearized equation but the wave equation for 
some axis lying (if possible) within the body. The 
approximations in pressure relation and angle of attack 
would again be avoided. 


Nonlinearity in lift has been shown to result pri- 
marily from viscous separation of the cross flow. The 
resulting force can be approximately accounted for 
using the simple viscous sweepback concept. Non- 
circular bodies could be treated similarly. A more de- 
tailed theoretical treatment of the cross-flow separation 


seems desirable. 
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Methods for Calculating the Flow in the 
Trefftz-Plane Behind Supersonic Wings’ 


P. A. LAGERSTROM?# ann MARTHA E. GRAHAM? 
Douglas Aircraft Company, Inc. 


SUMMARY 


By means of the linearized nonviscous theory, methods are 
developed, using the complex functions of conical flow, for find- 
ing downwash and sidewash in the Trefftz-plane—i.e., at an in- 
finite distance downstream, behind a supersonic wing. Equa- 
tions are derived which describe this flow field behind several 
types of ‘‘flat plate’’ wings: the wide delta wing, the trapezoidal 
wing (wide delta wing with cutoff tips), and the delta control 
surfaces (a narrow deltas wing, the two sides of which are de- 
flected differentially). Graphs are included which show the re- 


sults for these cases. 


As observed previously in connection with other cases, in the 
plane of the wing the flow field is characterized by infinite up- 
wash straight behind the wing tips and, in the case of the delta 
wings and control surfaces, by a region of constant downwash 
behind the wing. These conditions are modified immediately 
upon moving above or below the plane of the wing. It is found 
that the flow field behind the delta control surfaces, the two 
sides deflected equal amounts in the opposite direction, does not 
vary with changes in Mach Number. 


In applying the results presented here, one must, of course, 
take into consideration corrections due to effects of fuselage 
These effects are not, 


interference, nonlinearity, viscosity, etc. 


however, the subject matter of this paper. 
NOTATION 
Geometry of the Wing 


a = angle of attack 
= sweepback angle of leading edge 


¢ 
¢* = reduced sweepback angle of leading edge = arc tan f 
v = arccosf (see Fig. 3 or 8) 

m = Mach angle 

b = span 

b = (b/2)/cm 

c = root chord 

cr = chord at tip of trapezoidal wing 

f = mtang = tan ¢* 

m = tana = 1/V M? — 1 

A.R. aspect ratio 

A.R., = reduced aspect ratio = A.R.V M? — 1] 

x. taper ratio = tip chord/root chord 


For the trapezoidal wing, the parameters }, f, A.R.r, and T.R. 
are related by the following equations: 


A.R., r 4b ie ate 
b= ; (1+ T.R.) or A.R., ms where 1 S 6 S 
O<f <1 Ta 
Received June 12, 1950. 


* This paper was originally issued as a Douglas Aircraft Com- 
pany Report.! 

t Consultant, Santa Monica Plant. Also, 
fessor of Aeronautics, California Institute of Technology. 

t Aerodynamicist, Santa Monica Plant. 


Associate Pro- 


T.R. = 
A.R.r 2 


_ 4 1— Te.) 
~ AR, \L + TR. 


— of 
4/(1 + T.R.) 


Coordinates 
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= Cartesian coordinates of 
physical space with origin 
at intersection of leading 
edge and root chord of 
wing. x; positive toward 
right tip, x2 pesitive up- 
wards, x3 positive in free- 
stream direction 


x; VG = 1, 2, 3) 


‘ ( «a x2 | *) ‘ : 

t;( xX, = 2 = , ¥; = — }= reduced physical coordi- 

cm cm c nn 

R; = x;/(b6/2) = dimensionless coordinates 
used for narrow delta 
wing 

&; (7 = 1, 2) = *; or ¥; (for narrow delta 
wing) 

r,@0 = polar coordinates of plane 
x3 = c or Trefitz-plane 
(%3 = ©);7 a Vx? 2+ x? 
and @6 = arc tan Xo/a) 

? = r/om,? = r/(b/2) 

7 os". 9" = polar coordinates of Trefftz- 
plane with origin dis- 
placed along x,-axis [Sec- 
tion (3)]} 

a,v = polar coordinates of e-plane 

Va, 3/2 = polar coordinates of »v-plane 

7; (j = 1, 2, 3, 4), 06; (7 = 1,2) = notation in eplane or in 
y-plane (see Figs. 8 and 
14) 

s=i+t+ik= |z\e% (|3| =forr) = complex coordinate in 
Trefftz-plane 

« = ae” = complex coordinate in e- 
plane 

y= Va e9/2) = complex coordinate in p- 
plane 

Note: Primes are placed in the coordinates of the constant 


pressure wing to distinguish them from those of the wide delta 
wing when the two are being combined to give the trapezoidal 


wing [Section (4)]. 


Velocity 
M = free-stream Mach Number 
uj = perturbation velocity in x,-direction = sidewash 
us = perturbation velocity in x-direction = upwash 
u3; = perturbation velocity in x;-direction; is proportional to 
pressure 
U = complex sidewash function of e-plane 
Us = free-stream velocity 
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INTRODUCTION 


N A PREVIOUS REPORT,’ the induced flow field (side- 

wash and downwash) behind certain lifting wings 
in supersonic flow was studied with the aid of linearized 
nonviscous theory. Calculations were carried out for 
points at finite distances behind the wing and also for 
the so-called Trefftz-plane—i.e., a plane perpendicular 
to the free stream at an infinite distance downstream of 
the wing. By comparison it was found that results ob- 
tained for the Trefftz-plane could be used as a fair ap- 
proximation to the flow beyond a certain finite dis- 
tance behind the wing (reference 2, graphs 2 and 10). 

This is of considerable practical importance, since, in 
general, calculations for the Trefftz-plane are much 
simpler than those for points at finite distances. The 
following procedure is used for the Trefftz-plane. First, 
the sidewash at the trailing edges of the wing is found 
This sidewash remains constant in the free-stream direc- 
tion behind the wing and is thus known on the axis that 
in the Trefftz-plane corresponds to the plane of the 
wing. In the Trefftz-plane, sidewash and downwash 
may be combined to an analytic function of the com- 
plex coordinate in that plane. The problem is thus 
reduced to that of finding an analytic function whose 
real part coincides with the prescribed sidewash in the 
plane of the wing. From a mathematical point of view, 
this problem is exactly the same as the classical prob- 
lem of finding the flow around a two-dimensional lift- 
ing wing of prescribed lift distribution in incompressible 
flow. For certain problems the sidewash given is so 
simple that the potential problem may be solved by 
inspection. This was done in reference 2 for the flow be- 
hind a rectangular wing (page 59), a wing with raked 
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supersonic side edge and straight leading and trajj. 
ing edges (page 58), and the narrow delta wing (page 
60). 


For many important cases, however, the solution of 
the potential problem is no longer immediately obvious, 
The purpose of the present report is to develop a 
straightforward simple procedure for finding an ex. 
pression in closed form which describes the flow in the 
Trefftz-plane for a much larger class of wings than the 
above-mentioned types. It applies to wings of finite 
chord obtained from conical wings by cutting them by 
a line perpendicular to the free-stream direction (in 
other words, the trailing edge has no sweepback). Of 
course, this also entails the solution for wings that may 
be represented as the result of superposition of a finite 
number of cutoff conical wings. 


The method proposed in the present report consists 
of a direct substitution into the expression for the so- 
called complex sidewash function for conical flow (refer- 
ence 3, in particular, page 33). For most cases of prac- 
tical interest, this function is given explicitly in refer- 
ence 3 or may be easily computed from the data given 
there. 


For the purpose of illustrating the method and also 
because of its importance in itself, the wide delta wing 
is discussed in detail in Section (1). Results are given 
in Section (2) and in Figs. 4, 5, and 6. Another inter- 
esting case is that of “delta control surfaces’’—i.e., a 
narrow delta wing whose right and left halves are at 
opposite angles of attack. Here, the basic formula for 
the Trefftz-plane may be obtained by inspection, but 
the method discussed in Section (1) will prove valuable 
for the numerical evaluations. This is discussed in 
Section (3), and the results are shown in Fig. 10. Fi- 
nally, as an example of the superposition method, Sec- 
tion (4) contains a discussion of wings obtained from 
the wide delta wing by cutting off the side tips to ob- 
tain straight side edges. The methods developed in 
Section (1) will be applied to the components of this 
superposition. Results for a few particular cases are 
plotted in Figs. 5, 15, 16, and 17. 


(1) Wipe DELTA WING. DESCRIPTION OF METHOD 


Consider the flat lifting wide delta wing in Fig. 1. 
It will be convenient at first to assume that the wing 
root chord ¢ is unity and that M = 7/2. In the final 
formulas the results will be generalized to an arbitrary 
c and to any M for which the reduced aspect ratio 
AR. = AR.VM? — 1 > 4. 
A.R., means that the edges must remain supersonic. ) 


(This restriction on 


The sidewash 1; along the trailing edge at points right 
above the plane of the wing (x. = 0+) is given by the 
formulast 


+ These formulas may be derived by various methods. See, 
e.g., reference 3, Section II or the discussion of Eq. (5). 
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al, tan ¢ 
id er 
ei T — tan* ¢ 


x, — tang —x, — tan ¢ 
arc cos — arc cos (1) 
1 —x,tan¢ 1+ x, tan ¢ 


—al, tang 


1<x<ctng: wu = constant = % = mF 
1 — tan’ 9 


ane <2 — 1: he —h 
x, >|ctn g|: tm = 0 


For points right below the wing (x. = 0~), the sign 
of m is reversed by antisymmetry. Since the Mach 
waves from the trailing edge do not introduce any 
sidewash, Eq. (1) is valid both immediately upstream 
and immediately downstream of these Mach waves. 

The flow in the Trefftz-plane is then obtained by solv- 
ing the following potential problem: “, — 7 is an an- 
alytic function of the complex variable 2 = x; + 1x2 
satisfying the boundary conditions in Fig. 2. 

While there are standard methods for solving such a 
potential problem (e.g., a conformal mapping followed 
by use of Poisson’s integral), they do not give the 
answer readily in closed form. However, the Buse- 
mann method of conical flow suggests the following 
procedure. Consider the ¢-plane as introduced in ref- 
erence 4 (also discussed in reference 3). It is the 
plane x; = 1 (¢ = 1) with the complex coordinate 


edefined by 


id / =9 - 
e=ae, a=(1 — Vi —?f)/? 
, i = 7 ait 
7= Vi, + 2X7, & = arc tan (X2/%x1) (2 
2) 
I; = X%1/cm, Ze = Xo/cm 
Forc = 1, M = v7/2, 4: = 1, Ze = %. 


The Mach cone from the origin intersects the plane 
x; = 1 in a circle (‘‘Mach circle’). There exists an 
analytic function U(e) such that uw, = Re({U) inside 
the Mach circle. Usually, U is computed from the 
corresponding values for “2 or v3 (see examples in refer- 
ence 3). However, in this particular case it may be 
found directly from the boundary conditions in the e- 
plane (Fig. 3) by considering the singularities at the 
points ¢; (reference 3, page 60). 


The solution is thent 


—it,. (€ — &) (€ — &) 
U(e) = In (3) 
TT (€ —_ €3) (€ —_ €4) 
where 
i, = —al, tan g/V 1 — tan’ ¢ 


} These boundary conditions are actually the conditions for , 


the so-called symmetrical case. However, as discussed in refer- 
ence 3, page 74, the lifting case is deduced from this case by a 
procedure that essentially consists of using the same expression 
but by regarding OAe, as a slit and thus continuing the function 


from Q to Q’ as indicated by the dashed line. 








(mM=¥2) 
Q(©) =ae‘* / 








Boundary conditions in e-plane for U. 


Fic. 3. 


The function 4, — 7 satisfying the boundary con- 
ditions of Fig. 2 may now be constructed by the follow- 
ing simple procedure. With the aid of the inverted 


Joukowski transformation ft 


2 I bc OE nig 
=e¢e-+-ore = 
Zz € Zz 
D) ° 
e=ae", 2=i& + th (4) 
- x1 7 P 
& = 1% = yk = Lo = Xo/cm TT 
cm 


a one-to-one correspondence is set up between the in- 
terior of the Mach circle in the e-plane and the whole 
Then the value of u; — iu at a point z is 
equal—within an imaginary constant—to the value of 
the known function U(e) at the point ¢«, which corre- 
sponds to z according to Eq. (4). The imaginary con- 
stant to be added to the function thus obtained is chosen 
so that uw = O for z = o. Thus, the flow in the 
Trefftz-plane behind a flat lifting wide delta wing is 
given in parametric representation by the formula 


z-plane. 


. . o . y 
u; — iu, = analytic function of z 
—1t (e — &) (€ — €) | 
= E + | 
T (e — €3) (€ — &) | 
1+ V1 — tan’ g | | 
In 7 

1— V1 — tan’ (5) 
ad; = —al, tan o/V1 — tan’ ¢ 
2 1— Vi — 2 | 
=e+ -ore = a | 
Z € Zz 
Z2=>x+ 1X2 } 


t This transformation was used in reference 5 to obtain solu- 
Reference 3 uses (after reference 4) 
However, in 


tions in conical flow theory. 
analytical extension. This is a matter of taste. 
the present context, the use of this transformation appears na- 
tural. 

+t For the narrow delta wing, £ and & are defined in a slightly 


different way. See footnote, Section (3). 
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Fic. 4. Wide delta wing (f = 0.75): Velocity component normal 
to radii from origin in Trefftz-plane. 


A proof of the correctness of this procedure and a 
discussion of the evaluation of Eq. (5) will follow below. 
First, we wish to point out the following important 
difference between U as a function of € and m, — itz 
as a function of z. Both real and imaginary parts of 
u, — tz have physical meaning (as sidewash and down- 
wash) and this for every value of z. Furthermore, z 
is a physical coordinate. On the other hand, only the 
real part of U(e) has physical meaning and this only 
for ‘e| < 1 (where it means the sidewash inside and on 
the Mach cone in the plane x; = 1). Furthermore, e 
is not the physical coordinate in the plane x; = 1 but 
is connected with the physical coordinates x; and x2 
through the Tschaplygin transformation as given by 
Eq. (2). 

In proving Eq. (5), one should observe the following 
(a) For e real, 
Fur- 


properties of the transformation, Eq. (4): 
\e| < 1, 2is also real and < 1 in absolute value. 
thermore, in this case (but only in this case) the trans- 
formation « = (1 — vi = z*)/z and the Tschaplygin 
transformation, Eq. (2), coincide. (b) For e€ on the 
unit circle, z is also real and z > 1; for, if e = e’’, then 
z = sec J. (c) In particular, the point « = e” is 
mapped on z = ctn ¢, since cos y = tan ¢ (reference 3, 
page 42). 

Now the sidewash at the trailing edge for |x, < lis 
given by the real part of U(e) if physical coordinates are 
reintroduced according to Eq. (2). Since the sidewash 
is constant in the free-stream direction, it follows from 
(a) above that Eq. (5) satisfies the boundary conditions 
in the Trefftz-plane for |21| <A. 

For the region 1 < x; < ctn g, x2, = O*, Eq. (5) gives 
a value of “,; in the Trefftz-plane which is the same as 


-MARCH, 1951 
the real part of U(e) on the part of the Mach circ 
where 0 < 3 < y [remarks (b) and (c)]. 

On the Mach circle a = r—i.e., € represents the phys- 
ical coordinate in the plane x; = 1. Furthermore, as 
remarked above, the real part of U(e) represents the 
sidewash. But the sidewash on the part of the Mach 
circle 0 < 3 < ¥, x; = 1 is exactly the same as on the 
adjacent part of the wing outside the Mach cone 
namely, %;. Thus, the boundary conditions are satis- 
fied here. 

In the same way one shows that the boundary condi- 
tions are satisfied for the other points 21 o 1t == 

Finally, as is easily seen, the imaginary constant in 
Eq. (5) is adjusted so that uw. = Ofors = ~. 


(2) SUMMARY OF RESULTS FOR THE WIDE DELTA WING 


For the wide delta wing (Fig. 1) with sweepback = 


yg, chord length = c, free-stream velocity = Ul’. and 
free-stream Mach Number = WM, the sidewash x, and 


downwash #2 in a plane perpendicular to the free stream 
and at large distances behind the wing are given as 
functions of the coordinates x; and x2 by the formulas 


, auf ., (€ — a) (€ — &) 
Uy — WM. = 7 In = is 
V1 at i (€ — €3) (€ — «&) 
 1+Vi-f 
1 In - 7 (6a) 
l1-Vl-f 
€ = e”, €& = —¢*. ea = —e¥, 
g =e ™, y = arc cos f 


/ = , 
e = (1 — Vl — 2?)/2 or 2/2 = € + (1/e) (6b) 
z= & + 1% = (x1/cm) + 1(x2/cm) (6c) 


Here, the correction for chord length and Mach Num- 
ber has been achieved by multiplying x; and x2 by 
VM? —1/c=1 /(cm) and by using the reduced sweep- 
back angle g * defined by tan g * = tan ¢ tan yp = 
mtang =f. 

In evaluating these formulas, it should be remembered 
that m is antisymmetric with respect to both the x- 
axis and the x2-axis and “2 is symmetric with respect 
to the same axes. Thus, it is enough to consider the 
quadrant x; > 0, x2 > 0. 

In the plane of the wing Eqs. (6) reduce to the follow- 
ing formulas in nonparametric form: 


0<H#CI wm | f 
x =O0+taU. V1 —f 
( ti —f —F, =f) 
are cos —— =< @Fo coe a 
\ 1 — fX l + i vj - 
(4a) 
"19 4% < 1/f) / = 
ie = ors aU, = —f/V1-f 
mh > 1/f ; 
' thw al. = 0 ; 
x = QO" ) 
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METHODS 


CAS 1 wm - , 

0 S 71 t ae 7 x | 

XxX, = OT al. rv 1 —f 
l+vi1-f | 

In , 

l-vVi-f| 

z 21) me f | 
1 7 2 

pa a y, * 

= ots al, V1 —f (7b) 
V1 (1/#;7) + V 7 f2 | 
In : — | 
V1 — (1/#?) -V1l-f? | 
ae ae 
In : 
1-vi-f]] 

The value of “ for x, = O in the Trefftz-plane has 
been plotted in Fig. 4 for the value of the reduced 
sweepback angle such that m tan ¢ = 0.75 and in 
Fig. 5 for m tan g¢ = 0.5. The downwash is constant 


for points in the Trefftz-plane straight behind the part 
of the trailing edge which is inside the Mach cone. 
More precisely, in this interval, u2/alU., depends only 
on m tan g but not on x;. The variation with m tan ¢ 
is plotted in Fig. 6. Note that right at the trailing 
edge the downwash is variable inside the Mach cone 
but constant outside. This follows from the discussion 
in reference 2, pages 32 and 45 ff. 

In evaluating the flow for points not in the plane of 
the wing, it is convenient to keep the parametric form 
of Eqs. (6) and to use the following method both 
for numerical computation and for rapid qualitative 
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Wide delta wing (f = 0.5): Downwash along x,-axis of 
Trefftz-plane. 
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Fic. 6. Wide delta wing: Effect on downwash in Trefftz-plane 
(x1 < cm, x. = O) of varying reduced sweepback angle (*). 
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Fic. 7. Mapping from Trefftz-plane to auxiliary e-plane 
estimates. Given a point (P) in the z-plane, one 
first finds the corresponding point (Q) inside the 


This is described by Eq. 
Then it ts known that the 


unit-circle in the e-plane. 
(6b) and plotted in Fig. 7. 
real part of Ue) at (Q) is the sidewash at (P) and its 


negative imaginary part is the downwash. But L’(e) 
may be given a simple geometric interpretation. With 
the notation of Fig. 8S we have 
; al. f 
u,(Q) = Re{l (e)] : lo — 6) (Sa) 
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Fic. 9. Delta ‘‘control surfaces.” 
cos 0; (rn? + i — 4sin? p)/2nr, 
Cos 02 = (ro? + 73? — 4 sin? W)/2rers 


ux(Q) = —Im[U(e)] = 


ror 1+ v1 —f 
ie = f (Sb) 
rire l1— Vv ‘1 — f? 
r)° = 1 + a? — 2a cos (8 — p) 
ro” = 1 + a? + 2a cos (8 + py) 
r;?> = 1 + a® + 2a cos (8 — y) 
r4? 1 + a? — 2a cos (8 + yp) 
aand # are given in Fig. 7 or by the equations 
en - — OV rite sin [(¢1/2) + (¢, 2)] 9 
be + be? ” 
oe = are tan [&(1 + a?)/i(1 — a’)] 
2 £ 
y 4 Se 
¢i1 = are tan a ie es = 
& — | sin ¢ 
fy fo 
Se Se 
ge = arc tan ; ,=- Ps 
&+ 1 SIN ge 


al). f 


rV 1 — f? 


SCIENCES MARCH, 1951 
(for &, & in first quadrant, # is in first quadrant; fo, 
£1, & in second quadrant, # is in second quadrant, etc.), 
This geometrical interpretation may be used for an 
immediate rough estimate of the values and may also be 
used for more careful graphical or numerical computa. 
tion. Results for f m tan g = 0.75 have been 
plotted in Fig. 4.7 


APPLICATION TO OTHER CUTOFF CONICAL WINGs, 
DELTA CONTROL SURFACES 


(3) 


The methods described above may easily be applied 
to other cutoff conical wings. Some complications 
arise when the leading edge is swept upstream. How- 
ever, such wings will be used only as components in a 
superposition procedure, and, as will be discussed in the 
following section, the complications will disappear in 
the final result and may be avoided altogether by a suit- 
able choice of components. If the wings are cut along 
lines that are not perpendicular to the flow, 
serious modifications of the procedure will be necessary. 

As is obvious from the description of the method in 
Section (1), the flow in the Trefftz-plane may be evalu- 
ated as soon as the function U(e) is known for the cor- 
responding conical wing. This function has 
evaluated for a large class of wings in reference 3 and 
may easily be computed by the methods described and 
illustrated there. Some further examples will be 
treated in Section (4). Below will be discussed a case 
where the analytic function in the Trefftz-plane may be 
but the ideas developed 


more 


been 


obtained by inspection, 
Section (1) will prove convenient in evaluating the 
function. 

The example referred to just above is that of a narrow 
delta wing with a discontinuity in angle of attack at 
the center chord as shown in Fig. 9. This case would 
correspond to a pair of delta control surfaces deflected 
in opposite directions if the fuselage interference were 
neglected. 

From the material given in reference 2 it is easy to 
compute the sidewash at the trailing edge (see Appen- 


dix). The formulas are 
|, “=> - al ie x 
2 <= Viow 2 
n 
* £, V1 =n (10 


Eq. (10) is true for arbitrary chord and Mach Number. 
However, there is no dependence on these parameters} 
as long as the leading edges are subsonic. 

From Eq. (10) the analytic function in the Trefftz- 
plane may be derived by inspection by simply replacing 
Ral by z = %, + 7%. A constant is added so that as 


t This same graph is presented in more detail in reference 1. 

t It can easily be seen that this is consistent with the use of 
the Prandtl-Glauert transformation as described in reference 3, 
Section I.M. 
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zg © uy — 120. The result is 


My — Me = —al, x 
pees 2 
n i 7 
. s rV1 — 2 ‘ an) 


In evaluating the first term, a parametric representa- 
tion is used which corresponds exactly to the previously 
used inverted Joukowski transformation [2/2 = € + 
(1/e) ore = (1 — vWi- 2?)/s, where « = ae‘’].t 
Hence, this term may be evaluated by means of the 
transformation graphs (Fig. 7) or from Eq. (9). To 
evaluate the second term, it is convenient to introduce 
the following notation: 


b\? " Xe 
oe 2 * a 
= x1 — =) + Xe 6* = arc tan 
\( ; 2 7 i (b 2) 
| b 2 ’ 
ay v2 
re* = WZix, +.) +2”, 6** = arc tan 
Ve 2) s x1 + (b/2) 
p* = r*/(b/2) and #** = r**/(b/2) (cf. reference 2, 


page 61; cmty = 6/2). 
Then the real and imaginary parts of Eq. (11) give 
for the sidewash and downwash, respectively, 


2 > Q* g** 

u = —aU, E Ina + ~\/ papae sin 5 ) (12a) 
20 2 o* + o** 

tw = aU, (= —1+ a ppt cos 9 ) (12b) 


The following symmetry properties should be noted. 
Since the wing is assumed to be of zero thickness, 12 
is symmetric and u; antisymmetric with respect to the 
plane of the wing (x,x3;-plane). Because the boundary 
conditions in a are antisymmetric with respect to the 
vertical plane through the center chord (x2x;-plane), it 
can easily be seen that uw is antisymmetric and 
symmetric with respect to that plane. Hence, it is 
only necessary to evaluate Eqs. (12) for the first quad- 
rant (x;,X%2 2 0). 

By means of Eqs. (12) and (13) the velocity com- 
ponent normal to radii from the origin of the Trefftz- 
plane is plotted in Fig. 10. 

In the plane of the wing, Eqs. (12) reduce to the fol- 
lowing formulas in nonparametric form: 


0< x <, B Xo = Ot, uy, = —al, 4 
2 1-V1—#? 2 
- In z + F, = :) (13a) 
1 rv 1 rl” 


( 2 l 2 
— — are sin — + 
v x rv &2 ac +) 


Tt Eq. (4), where now s = & + t& = ¥ + 7% = x,/(b/2) + 
i[x2/(b/2)] 
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Fic. 10. Delta ‘‘control surfaces’: Velocity component normal 


to radii from origin in Trefftz-plane. 
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Fic. 11. Trapezoidal wing. 


These formulas may then be completed for negative 
values of #; and %, according to the symmetry laws de- 
scribed above. Notice that, as in practically all ex- 
amples studied, the curve of uw vs. %; for % = 0 hasa 
flat part immediately behind the wing. The present case 
is distinguished by the fact that there is one flat part 
for 0 < &, < 1 and then (because of antisymmetry) 
a discontinuous jump to another flat part for —1 < 
x, < 0; furthermore, the value assumed there is inde- 


pendent of Mach Number. 


(4) SUPERPOSITION METHOD. TRAPEZOIDAL WING 


It is obvious that by the method developed in the 
preceding sections one may obtain the flow in the 
Trefftz-plane behind wings that may be represented 
as a superposition of cutoff conical wings. As an ex- 
ample the trapezoidal wing in Fig. 11 will be discussed. 
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Fic. 12a (top). Trapezoidal wing tip. 
Fic. 12b (bottom). Half infinite sweptback wing. 
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Fic. 13. Constant pressure wing adjacent to flat plate. 














¥ 
va) 
Re 
e 
/3 
eg 
% Z 
Fic. 14. v-plane: Notation for Eq. (15). 


There are many ways of expressing the superposition. 
The basic components for a simple representation will 
be related to the wide delta and the wings in Fig. 12. 
Actually the conical wing obtained by subtracting the 
wing in Fig. 12a from the wing in Fig. 12b will be used. 
It will consist of a triangular part DEB of constant 
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Fic. 16. Auxiliary constant pressure wing (f = 0.5): Downwash 
along x-axis of Trefftz-plane. 


pressure (but variable a) adjacent to a flat plate at a = 
0 (see Fig. 13). The constant pressure wing induces a 
lift distribution on the flat plate inside the Mach cone 
from D. 

The trapezoidal wing (Fig. 11) is now represented as 
a wide delta wing whose tips are cut off by subtracting 
wings of the type shown in Fig. 13. 

The solutions in Figs. 12a and 13 differ only by a con 
stant (obtained from the wing in Fig. 12b). However, 
the wing in Fig. 13 has the advantage of being effec- 
tively of finite span (all perturbations vanish to the left 
of the Mach cone from D). For this wing, the con- 
struction in Section (1) can be carried through without 
any difficulties. For the wing in Fig. 12a, the same 
construction would yield sidewash and a singularity 
even to the right of the point in the Trefftz-plane which 
is straight behind -. This is obviously an incorrect 
solution, and, hence, it is undesirable to use it even if 
it turns out that after the complete superposition has 
been carried out the errors compensate.f Also the 
algebraic subtraction suggested corresponds in an in- 
tuitive way to the notion of actually cutting off the 
wing tips of the wide delta wing. 

+ Asa matter of fact, the solution for the wing in Fig. 12b, which 
also has an infinite span, would exhibit similar undesirable fea- 
tures. If the solution for the wing in Fig. 12b is subtracted from 
the one in Fig. 12a, the errors are subtracted out and the correct 
solution for the wing in Fig. 13 is obtained. However, it is 
easier and safer to work with the conical wing in Fig. 13 directly. 
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The complex sidewash function is obtained by subtracting from a constant Eq. (3.9) in reference 3. It is 
- +241 | y ils (vy — wm) (Pv — M%) 
U(e) = y+ cos — + cos wy In + iy (14) 
T v 2 T (v — vo) (v — v3) 
plus an imaginary constant, where 
— ol i f al & 


_ ’ 113 / 
V1i—-f Vi-f 


As was the case for the wide delta wing, the parametric representation above suggests the following computa- 
tional procedure. A mapping will be constructed which takes a point P = x, + ix. in the Trefftz-plane into a 
point Q in the complex v-plane. Formulas will be given for evaluating ™ and u2 as functions of Q in such a way 
that their values at Q represent the true value at P. The mapping is defined as follows: 

The point P = x; + ix: is mapped on the point ¢ = ae’’ in the e-plane by the inverted Joukowski transformation, 
Eq. (4). Q is then defined as the point ~/ae’’”. Stripped of complex notation, the above definition states that 


if P has Cartesian coordinates x), X2 then Q has polar coordinates ~/a, 3/2, where a and # are defined by Eq. (9) 


or by Fig. 7. The formulas for m4; and 1, in the v-plane are, then, 


wea : V2 , v l f 
u(Q) = Re[U(e)] = —al, al _, (sin 5 (. a ne + alt _ (0. — +n) | (15a) 


: : Vv ‘2 / v ; l I riv4 
u(Q) = —Im[U(e)] = —aU, V1 _, (os Sv a+ = + eer Sard pr In “i - 
2 f Vitf~1 ’ 
] — =In 5 = | (15b) 
mVi-f rVi-fi Vitft 


6; and 7; are shown in Fig. 14 and may be evaluated in the same manner as described for the wide delta wing. 


Analytically, they are defined as follows: 


6 = arecos {[r? + re? — 4 cos? (/2)]/2nre} 
6. = are cos {[r3? + ry — 4 cos? (W/2) |/2r374} 
rn? =1+a— 2vV/acos [(W/2) — (8/2)] 
re? = 1+a+ 2v/acos [(y/2) + (8/2)] 
r??> = 1+a+t 2Vacos [(~/2) — (8/2)] 
rn? = 1+ a — 2vV/acos [(y/2) + (8/2)] 


In the plane of the wing one may express “4 and wu directly as functions of x; and A». 


—~1¢ 4 <0) j__2 l+ a f ene teN | 
+m = aUy os - / .. arc COs 
Xe = Ot Urv1 -f — i; rVvil—f? 1 — fz, f | 
0c aC 1s : 
oe = po uw, = —al,, (f/V 1 — f?) (16a) 
a | 
A<-L&> I/fl, _g 
ve = Oth 
-1<% <0) 2 Vi+f } 
, +¢ 2 = al, + = In | 
x2 = 0 eV1i-f 2eV1-f V14+74+1 | 
i<—1%>0 ff 8 / + (16 
ao of me = aU 1- V+ &)/m] - f <oeey 
x. = 0 Invi —f 
f Vi+tf—-WVW(1+ #)/a 4 Vi4s~ 18 | 
7 n 7 sete ie > | 
mrVil-f V1i+f+V(1+ &)/t Vitf+i1)) ) 


where %; = x,/cm, c = DE of Fig. 13. 
The downwash in the plane of the wing for the particular constant pressure wing with f 
The variation of downwash with reduced sweepback angles for the region —1 


= mtan g = 0.5 is 
€£ % < Ois 


plotted in Fig. 15. 
also shown on this figure. 
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The downwash and sidewash for the trapezoidal wing is then that given by Eqs. (7) or (8) for the wide delta wing 
minus that given by Eqs. (15) or (16) for the wings of constant pressure adjacent to a flat plate. For any trape. 
zoidal wing for which m tan g = 0.5, the downwash in the plane of the wing may be obtained by combining the 
plots of Figs. 5 and 15. Each of these graphs involves a logarithmic infinity at the point corresponding to the out- 
board tip of the original delta wing, but this disappears when the two are combined. In order to facilitate this com- 
bination, Eqs. (7b) and (16b) have been modified so that from each this singularity is eliminated. This is done 
by removing the logarithmic term from Eq. (16b) and by subtracting this same term from Eq. (7b). While 
neither has any particular physical significance alone, their difference still gives the downwash for the trapezoidal 
wing, since both formulas are modified by the same amount. 

In combining these formulas, <, of Eq. (16b) will be primed to distinguish it from %, of Eq. (7b), since they are 
now related. Single primes distinguish coordinates of the right tip constant pressure wing; double primes, the 
left tip constant pressure wing (see below). 

Eq. (7b) modified is 


oe SL ivi -f +Vi-@ |, |\Vits— Vi + 1 | 
aU. wrV1i-fi| aie * "Vit f+ Vi + [y (#,')]| 


Vits—-Vi+ Wei _ | LEVIES) — ag |t= VEEL 
IVi tft V1 + [1/(%]| j~—-Vi-f] + vi¢f4i\ 


= gi(%1) + go(%1’) + gs(E1") + constants (17) 


In 


Note: for —1 < % < + 1 omit gi(%:); for —1 < %’ < 0 omit go(%’); for —1 < %” < 0 omit g;(#,") 


J i xy" x, — (b/2) %i-b _, x3" —x, — (b/2) —# —b 
ti = xi/em, i,’ = = = - ff a—— = ~ a 
Crm cm T.R. 1 — Of Crm cm T.R, 1 — Of 
where c = root chord, 6 = span, and T.R. = taper ratio of trapezoidal wing. cr = tip chord of trapezoidal wing 


= root chord of constant pressure wing. Eq. (16b) for the constant pressure wing to be removed from the right 
tip is, when modified, 


ue/aU. = (2 nV -fitl- v1 + (1/z')] (18) 


Eq. (18) holds for the left tip if 7,’ is replaced by %,” as defined above. 

Eq. (18) is plotted in Fig. 15, and Eq. (17) is plotted in Fig. 5 for a few particular values of the parameter } = 
(b/2)/cm. It is to be noted that, for 6 small, corresponding to low reduced aspect ratios, the contribution to the 
downwash by Eq. (17) is small. 

The downwash at the exact point at which the singularity originally occurred, #; = 1/f, can be found from Eq. 
(17) only by a limiting process, though the value at this point may be found readily by interpolation. At # = 
1/f, Eq. (17) reduces to 


Beas tn)... Wide viz 
= {in a n / 
ae #V1i-f{ [Al —- #) VitftV1+ (1/H”) 
li-Vi+f t+ VIF} 
21n — 1n|- * (19) 
1+ Vit+f add tad i 


where %,” = (1 + Of)/—f(1 — bf); 1 <b < 1/. 


Again it is sufficient to consider only positive values Figs. 16 and 17. In Fig. 16 the curves for the rec- 
of Z, since “2 is symmetric about #, = 0. tangular and the trapezoidal wings show a marked 

The superposition process described above has _ similarity. At the same a, MW, and U., the lift of the 
actually been carried out for two particular trapezoidal rectangular wing is 0.677 times that of the trapezoidal 
wings, one with m tan ¢ = 0.5, h = 1, which corresponds wing (this factor is due to an area ratio of 2/3 and an in- 
to a reduced aspect ratio of 2.67 and taper ratio of 1/2 significant difference in lift coefficient). If the curve 
and the other with m tan g = 0.5, 6 = 1.5, correspond- for the trapezoidal wing is multiplied by this factor, 
ing to A.R., = 4.80 and T.R. = 1/4. The results are the two curves nearly coincide. However, the trape- 
shown in Figs. 16 and 17. zoidal wing shows a little more downwash at the center 

For comparison, the downwash of a rectangular wing and a little less at the tip. This is due to the fact that 
with the same span and tip chord has been plotted on for the rectangular wing there is no sidewash behind 
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the center section, so that the downwash around x; = 0 
js due only to the vortex sheet (sidewash discontinuity ) 
behind the tip Mach cone. The trapezoidal wing, on 
the other hand, has sidewash behind the center section, 
and this effect is more pronounced the larger ¢ is. 

Fig. 17 shows a greater discrepancy between the two 
curves. This is due to the geometry of the wings 
which can be seen from the plan forms that are drawn 
to scale for M = +/2. The rectangular wing has a re- 
duced aspect ratio of 12. Hence, the vorticity is con- 
centrated in a narrow region behind the wing tips. 
For the trapezoidal wing, on the other hand, the side 
edge effect is small. The downwash distribution is 
rather similar to that of the wide delta wing (Fig. 5). 
For the wide delta wing the value of mw. starts to in- 
crease to infinity when #; increases from | to its value 
at the tip. The beginning of this rise may be seen in 
Fig. 17; however, because of the influence of the side 
edge, the increase is checked for points behind the tip 
Mach cone. 

The irregular behavior of the downwash curve near 
the tip as shown in Fig. 17 is typical for aspect ratios 
for which the theory is valid. As A.R., decreases to its 
limiting value, it disappears gradually and is entirely 
absent in the limiting case (A.R., = 2.67 for m tan 
= (0.5). Itis this limiting case that is plotted in Fig. 
16. 


Appendix 


SIDEWASH AT TRAILING EDGE OF DELTA CONTROL 


SURFACE 
In reference 3 the complex pressure function W(e) was 
derived for the wing of Fig. 18. 


3) (B 
We) = ry + ie + 5) (A-1) 


(e — B) (Be — 1) 

where ¢ is defined in Eq. (2), B = (1 — V1 — b?)/b, 
b = (b/2)/cm, and k is a constant dependent on a; and 
a. This formula is valid only if a, and ay satisfy a 
certain relation.{ For e real, Re [W(e)] = us which is 
proportional to the pressure distribution on the wing. 

The wing of Fig. 18 has a discontinuity in slope along 
the x;-axis given by [see reference 3, Eq. (1.42c) ] 


— Aus —T (>) 
= = = > ; = 
‘Ie Oe a 2 ae 


T —kr 


_— l 
1(0+3)-, 
2U.m ( r B bU m 


Substituting Eq. (A-2) in Eq. (A-1) 


Wie) = mbU (a2 — a) (e + B) (Be + 1) (A-3) 
7 (e — B) (Be — 1) 


(A-2) 


Tt See footnote on page 190. 
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In a similar manner, for a wing such as that of Fig. 
18, but with a; and a» reversed, 

mbU..(a2 — a) € ae B i. 1 

ie ( 1 (e — B) (Be — 1) 

7 (e + B) (Be + 1) 

But the wing in which we are interested is that of Fig. 9. 

It may be obtained simply by letting ag — a, = a@ and 

subtracting from the wing of Fig. 18 the wing, with az 

and a; reversed.f This gives the angle-of-attack dis- 

tribution of Fig. 9, for which the complex pressure 


(A-4) 


function is 


mabU (e + B) (Be + 1) 7 
(e — B) (Be — 1) 


(e — B) (Be — 1) 
(e + B) (Be + 1) 
mabU ., (1 + B?) (2e) 


aa (A-5) 
s | VW(e — B?) (B%? — 1) 


We) = 


us 


t If the two wings had been added instead, the a; and a, couid 
have been determined as functions of k from Eq. (A-2) and refer- 
ence 3, Section III.J. However, these functions are of no interest 
in the present context. 
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Then, on the wing [e real = +a where 2/|%| = a+ 
(1/a)] at the trailing edge, the pressure perturbation js 
2mabU .¥, 


Re[W(e)] = —_————= 
aV 5? — #2} (A4) 


for || < b: U3 
for |z,| >b: u=0 


But here the following relation exists between the side- 


wash and pressure perturbation [reference 3, Eq, 
(1.43a) ]: 
du;/di, = (—1/mz,) (du;/dz;) (A-7) 
Hence, 
—2ab?U,, —2aU, 
du; - — ey di, = es > parse dx, (A-8) 
1X, (b? — dy) *” wX;(1 —_ X 1") s 
where 
xy = XxX /(b 2) = X ib 
or the sidewash is (within a constant) 
—2al, 
for [%,| <1: 4 = xX 
T 
~ 1 — a8 l 
he l l i, =~ (10) 
|%,| Vi1i-—- #,* 


for |*,| ~ i: ee = 
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Purely Rolling Oscillations of a Rectangular 
Wing in Supersonic Flow 


TING-YI LI* 
California Institute of Technology 


SUMMARY 


The purpose of the present paper is to show a method of solu- 
tion of the purely rolling oscillations of a rectangular wing in 
supersonic flow and to present some numerical results. From 
these results it is possible to construct vector diagrams of the 
rolling moments acting on wings of several aspect ratios (A.R. = 
1/8, 1, 3, 5, and ©) at two Mach Numbers (M = 10/7, and 2) for 
a range of reduced frequencies (0 < k, < 2.0). Thence, it is 
found that the tips of a finite wing tend to produce two effects— 
namely, (1) to reduce thé real part of the rolling moment vectors 
and (2) to shift the rolling moment vectors in the positive direc- 
tion. It is believed that these effects are of great importance for 
estimations of the damping in roll. 


INTRODUCTION 


- REFERENCE |, PRoF. H. J. STEWART and the present 
author developed a method of analysis for the bound- 
ary value problems of the vertical and pitching oscilla- 
tions of a rectangular wing at a supersonic speed. It 
was shown that the effects of the wing tips are to shift 
the phase of the lift and pitching moment vectors in 
the positive (leading) direction. In some cases, the 
phase angles of these vectors were found to be positive 
i.e., these vectors being leading the corresponding quasi- 
steady vectors. It is of interest to carry on further 
this kind of analysis to another mode of oscillation 
namely, the rolling oscillation. The purpose of the 
present note is to show a method of solution of the 
problem of purely rolling oscillations of a rectangular 
wing at a supersonic speed and to present some numeri- 
cal results. 

A flat rectangular wing of zero thickness placed in a 
main stream of uniform supersonic speed L/’ is supposed 
to perform a slight periodic motion such that on the 
top wing surface (Fig. 1) 


r+ ( 


l "No - : ivl / . 
d = exp (zvt) exp | — —, (x — &) | dé 
T 0 B-a- oy =e 


n 


where the area of integration for the surface integral is 
evidently that portion of the plane z = 0 intercepted 
by the forward Mach cone extended from the point P 
and the source strength has been determined from the 
given downwash condition in Eq. (1). 
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lim (O@/0z) = —UAgy exp (trt) (1)t 
z=+0 
It is clear that this simply represents a slight periodic 
rolling motion about the x-axis. The problem is to find 
the rolling moment about the x-axis, acting upon the 
wing surface at every instant, due to this oscillatory 
motion. 

The general procedure in reference 1 will be followed 
here. 

Firstly, the velocity potential functions in the purely 
and mixed supersonic regions will be separately deter- 
Secondly, the rolling moments due to the 
These 


mined. 
various regions of the wing will be computed. 
components will be summed up to give the total rolling 
moment. In the final presentation, it is convenient to 
use the dimensionless rolling moment coefficient, C,, 


defined as 
C, = R/(1/2)pU*b*x, (2) 


where R represents the total rolling moment about the 
x-axis. 
SUPERSONIC 


POTENTIAL IN THE PURELY 


REGION 


VELOCITY 


It is shown in reference 2 that for a periodically os- 
cillating wing at supersonic speeds it is possible to re- 
place the actual wing by a distribution of oscillating 
sources in the plane z = 0. Thus, the velocity poten- 
tial at a point on the wing can be expressed as a surface 
integral. For a point P(x, y, +0) in the purely super- 
sonic region on the top surface of a rectangular wing 
with purely rolling oscillations (Fig. 1), the velocity 


potential is 


VA ~ 9 9 cos (v/B7a) [(x — &)? — By — n)*]" dn 


aie - 8 (@ — 8? — By — a)” 


By introducing a set of new integration variables de- 
fined as 
a = (vU//B’a*) (x — &) | 
» (4) 
o cos 6 = (vU/Ba*) (vy — n) | 
+t Most symbols used in the present note are defined in refer- 


ence 1. 
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the integration with respect to 6 can be carried out, and 


Eq. (3) becomes 


Apa’ 


Py = : exp (ivt)yT(k, M) (5) 
V 


where, as in reference 1, 
k 
T(k, M) = f exp(—t0)Jo(a/M)do (6) 
0 


and 
k = (vx/U) (M?/B?) (reduced frequency) (7) 


d, = 


By converting into the variables o and 8, it is found that 


AERONAUTICAL 


UA x + ply — (6/2)] iyU 
 exp(ivt) f exp | - oe -e- | xX 
7 J0 B*a* 


f + (1/B8)(x — &) 
b — y — (1/B)(x— &) 


SCIENCES—MARCH, 1951 


IN THE MIXED SUPERSONIC 


REGION 


VELOCITY POTENTIAL 


In reference 2, a simple ‘‘equivalent area’’ theorem 
has been established to determine the appropriate area 
of integration for the surface integral representing the 
velocity potential at a point in the mixed supersonic 
region on a finite oscillating wing. By this theorem jt 
is known that the shaded area in Fig. 2 will not contrib- 
ute anything to the velocity potential at the point 
Q(x, ¥y, +0) on a rectangular wing. Therefore, the ve- 
locity potential at the point Q may be written as 


P; = Dy — Py (8) 


where 


n cos (v/B?a)[(x — £)? — By — n)?I 


7, — dyn (9) 
[(x — 8)? — 6%(y — n)*]” : 


UAy x ; , 7 . xo o 
d = - exp(zvt) exp(—ta)da y- cos 8 } cos sin@}d6 (10 
7B k -/ k (B/x)[(b/2) — y] cos~!{1 — [28/x(o/k)][(b/2) — y]} B k M 


This may also be written as 
eT 


d, = 


where 


= 


ROLLING MOMENT COEFFICIENT 
The rolling moments about the x-axis due to the 
three regions So, S;, and S2 will be designated as Ro, Ri, 
and Ro», respectively (Fig. 3). Then, it is easily seen 
that 


(6/2) — (x,/B) Ee Ody 
Ry = 2p f ydy f (ist + U ax (13) 
0 J0 Ox 


rb/2 


R, => 2p / y dy p 4 
(b/2) — (x,/B) 


¥e ‘ , OP; 
tvP, + Ll Ox (14) 

pl (b 2) — y] Ox 

b/2 
R, = 2p f ydy X 
(6/2) — (x ./B) 

‘i 2) — y] oO, 
ivy + L’ ) Ix (15) 
J 0 ( : Ox 


where %) and %, are as defined in Eqs. (5) and (8), 
respectively. 

The total rolling moment about the x-axis due to the 
complete rectangular wing is consequently expressible 


UA , : : . x k 
ra exp (tvt) y exp(—tky)du / (» ws ine 0) cos (= sin a)ao (11) 
7B (8/x)[(b/2) — y] e/ cos! {1 — (28/xp)[(b/2) — y]}} B M 


a/k (12) 


as 


R = 2(Ry + Ri + R2) 


6/2 *. : OP, 
4p y dy ivPy + U dx — 
Jo 0 Ox 
b/2 
4p f y dy X 
(6/2) — (x, B) 
Xe OP, 
[ (iva, +U ‘) dx (16 
pt (6/2) — y] Ox 


where % is as defined in Eg. (11). Finally, by 
Eq. (2), the rolling moment coefficient C, may be written 


as 
C, = C, — AC, (17 


where, from Eq. (16), 


8 b/2 
C, = ) ly 
Dx, ri dla 
"#e Ov; r 
/ (iv + U ) dx (18) 
J/ 0 Ox 


and 
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Fic. 1. A rectangular wing in purely rolling oscillations. 


b/2 
AC, = : J y dy X o, 
"Ube. S62) - (x,/B) 4 . 
>» © 


¥e ; _ OP» 
f (iva, + U ) dx (19) 
Bl (b/2) — ¥] Ox f 


i (x, ) 











Of course, C,, is the rolling moment coefficient of the 
rectangular wing computed on the basis of the velocity 





potential in the purely supersonic region, and AC, is , 
the necessary correction that accounts for the tip effects x 
of the rectangular wing. Fic. 3. The right half of the rectangular wing. 


EVALUATION OF C,, 


Combining Eqs. (5) and (18) yields the following result: 


2 . *b/2 %, b/2 
C,, = so exp(ivt) ot I y? dy f T(k,M)dx + T(k,,M) f i dyt (20) 
vUb?x, ‘Je ° ° 0 


The y-integration indicated here may be immediately carried out, and the variable x in the formula may be elimi- 


nated by the reduced frequency k by Eq. (7). Thence, it is obtained that 


C,, = (bAo/3B8k-) exp(ivt) [T(k.,M) + 1(8?/M*)A(k.,M) | (21) 


k 
A(k,,M) = [ T(k,M) dk 
/0 


as defined in reference 1. Introduce the expression [cf. Eq. (26), reference 1] 


where 


(1/2)bC,,.* = (2bA0/B) exp(trt) (22) 


which may be interpreted as the two-dimensional quasi-steady lift coefficient based on the instantaneous angle of 


attack at the right wing tip. Then, it is easily found that 
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Cr, 
12 C,,* 


B? 
M 


C, l ; 
—_ = | Tbe) + $ 


5 
bC,,* 12k; (23) 


| Ako) ~ 


where C,,/C,,* has been computed in reference 1 for a range of reduced frequencies (0 < k, < 2.0) at two Mach 
Numbers (7 = 10/7 and 2). Asinreference 1, Eq. (23) holds only for BA.R. 2 1. 
EVALUATION OF AC, 


Eq. (19) may be rewritten as 


8 iv *e a - 
AC, = = - | 4 f ax f yPo(x,y)dy + J yPo(x,,y)dy (24) 
l bx, l 0 /2) — (x/) (6/2) — (x,/B) 


where #,(x, y) is the function defined in Eq. (13). It is convenient first to deal with the integral 


b/2 —— b/2 1 
l Aox . P 
yPo(x,y)dy = — exp(tvi) y dy exp(—tku)du X 
b 7B (b/2) — (x/B) (B/x)[(b/2) — y] 


/2) — (x/B) 


ri x k UAox ; */B /b 
J (» ee s) cos ( * sin ) dé = exp(tvt) Fd ( _ y' Jay’ x 
cos~!{1 — (28/xp)[(b/2) — y]} B M 7B 0 2 
' 4 b x k 
f exp( ~ituddn f ( —-y— cos s) cos (; sin s) d@ (25) 
(B/x)y’ cos~![1 — (28/xp)y’] \“ B M 


where y’ = (b/2) — y._ In Eq. (25), by twice inverting the order of integrations, it is obtained that 


b/2 UA wt if ku, 
yPo(x,y)dy = exp(ivt) exp(—tku)du cos sin 0 ) dé X 
(b/2) — (x/8) TB 0 0 M 
*(xp/2B)(1 — cos @) b b x 
| ( _ y') ( —y'- e cos 0) dy’ (26) 
° 2 2 8 


Performing the integration with respect to y’ yields 


(xp/2B)(1 — cos @) b b x x bx 22 
( = y’) ( -y'- : cos 0) dy’ = : 6? — - tu ) + 
0 2 7 2 B 86 28 6B? 
pre? bxu = x*p? Pe 
—bh? + ) cos 6 ( _ ) cos 26 cos 30 (27) 
( opt ) 8° F \ 5g — age ) 98 % T Gaz °° 


After substitution of Eq. (27) into Eq. (26), the 6-integration may be carried out. The final result may be written 


as 
b/2 y e x 
UAc x ;  £ bse". . 
G(k,M) = yP,(x,y)\dy = ——- exp(trt) | b? - Lo (k,M) — L,?(k,M) — 
(b/2) — (x/B) 8B? k k 28 \k 
: G) 2 (k,M) + 2 ()1 *(k,M) — — (:)z 8(k,M | (28 
6a? \k/ 2a\k) * ©” ower r” 
where 


ek 
Lia” (k,M) = J o” J, (¢/M) exp(—ia)do (29 
0 


Introduction of the G-function in Eq. (24) yields 


+k 


AC, = (8/Ub*x,) [G(k.,.M) + 7(8?/M?) / 'G(k,M)dk] (30) 


0 


- 


It is interesting to notice the formal parallelism between Eqs. (21) and (30), but for the purpose of computations, it 


is desirable to use the following formula instead: 








di 


bo 


(26) 


ritten 


(28) 


(29) 


(30) 


1S, it 
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9 wk 
acre 1s 1 gE . a c | 1 
—_ = < ) k.,M) pee Lo'(k,M) 1k — 
bC,,* tke (ke BA.R. L™ oe J 2(keBA.R.)? * 


9 


8 Re 1 : ef, | 1 
2 ; Lo%(k, M)dk | — —aie S| : Ly3(k, M)dk 
[ ell) + ton I or | 6(k-BA.R.)? | enti’ I re 36a)? * 


ee ff”. | 1 [ . 8 f | | 
Ls*(k,,M) L,*(k,M)dk | — L.3(k.,M L.*(k,M)dk {+ (31) 
| ms + tap f (Me 2(k-BA.R.)3 oto. eer 


Eq. (31) is valid for BA.R. 2 1. 
REDUCTION FORMULAS 


From the analysis of reference 1, it is immediately seen that 


Lo'(k,M) = RT(k,M) — A(k,M) (32) 
Lo(k,M) = k?T(k,M) — 2B(k,M) (33) 
Ly'(k,M) = R®T(kR,M) — 3E(k,M) (34) 
k. 
‘ [ Lo'(k,M)dk = 2B(k,,M) — k,A(k,,M) (35) 
J 0 
*k. 
I Lyo?(k,M)dk = 3E(k,,M) — 2k,B(k,,M) (36) 
0 
k. 
f L,3(k, M)dk = 4F(k,,M) — 3k,E(k,,M) (37) 
0 


where, as in reference 1, 
k, k, "ke 
B(k.,M) = f kRT(k,M)dk, E(k,,M) = f k?T(k,M)dk, F(k,,M) = | k®7(k,M)dk 
J0 0 0 


Another group of reduction formulas are derived by integrations by parts. They are as follows: 


k l k l 
L2?(k,M) = exp(—tk) | tk*J2 | — *—1ak) J — — [Lo(k,M) — iL,'(k,M)} 38 
2?( exp(—ik) [ 9 (*) + uw tk) Ji (j,) TE | iL, ) | (38) 


k I k 
Li(k,M) me exp(—ik) (ie + k?) Jo () — Vl (—k’ + 31k? + 3k) J; (*,) | + 


l 
ps [— LotlkeM) + BiLo%(k,M) + 3Lo'(k,M)] (39) 


k 
a 1 3 : Baer 
[ L.*(k,M)dk = exp (—ik,) | (ik,2 + 2k.) J; ( )| 4 iLs%(ke,M) — — iLg%(ke,M) + 2L9'(keM)] — 
, M M WG 
1 -” i k. 
L,2(k,M)dk + ~- | ‘Lo(k,M)dk (40 
ra ee f mann 


k 
| ke ae 
[ L;*(k,M)dk = exp(—ik-) F , (ike + Ske — Sik. F )| + iln*(k,M) + Lo%(kyM) — 
0 d d 


I I Ms 3 Of* 
- [iLo?(RkepM) + 5L02(k.,M) — SiLo'(ke,M)] — — f Ly?(k,M)dk + 4 f Lo?(k,M)dk + 
M? M? Jo M? Jo 


3 we 
| f Lo'(k,M)dk (41) 
M? Jo 


Therefore, by these reduction formulas, it is possible to express the important functions that appear in Eq. (31), 
directly or indirectly, in terms of the functions 7, A, B, E, and F, which have already been computed in the prep- 


va 
< 


aration of reference 1, for a range of reduced frequencies (0 < k, < 2.0) at two Mach Numbers (M/ = 10/7 and 
2.0) 
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TAYLOR’S SERIES EXPANSIONS 


As in reference 1, it is interesting to find the first two terms of the Taylor’s series for the Eqs. (23) and (31), 
respectively, as follows: 


C,,/bC,,* = (1/12) [1 — (ék,/2M2) + O(R-2)] (49 
“a. 1 ([! Bi ll ]+ie bo toy 1 gi. 
bC,* 4BA.R.\|2  6BA.R. 24(BA.R.)?] 3 -8BA.R. | 30(BA.R.)2 M2 


] l 
E _ —- |} + o(k-)) (43 
4B8A.R. 20(BA.R.)? 
Combining Eqs. (42) and (43) yields 


CG, _C,—AC,_ 1 | 3 
<* «,* 


l 
. 1+ 
2BA.R. * 2(BA.R.)? ° aE 

tk M*+1 2M? + 1 3M? + 1 
“_ | - 1 — _— - 1 + OCR.) (44 

24M? BA.R. 4(BA.R.)? 20(BA.R.)* 
From these equations, it can easily be seen that, for small reduced frequencies, the real component of the rolling 
moment vector acting on a rectangular wing with finite span is usually considerably smaller than the real com- 


ponent of the corresponding vector acting on a similar wing with infinite span—for instance, when A.R. = 1/8, 
(C,/bC,z,*), = (1, 8) (C../8C;."), 


Furthermore, for the same small reduced frequencies, it is evident that, from Eq. (42), for an infinite wing the roll- 
ing moment vector will have negative (lagging) phase angle and that, on the other hand, from Eq. (44) for a finite 
wing the corresponding vector may swing forward to have a positive (leading) phase angle at certain appropriate 
combinations of Mach Numbers and aspect ratios. Therefore, the effects of the wing tips are (1) to reduce the 
real part of the rolling moment vectors and (2) to shift the rolling moment vectors in the positive direction. These 





effects, of course, can greatly influence the calculations of the damping in roll. 
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Fic. 4. Effects of the aspect ratio on the rolling moment vector, C,/bC1,*, at M = 10/7 
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Fic. 5. Effects of the aspect ratio on the rolling moment vector, C,/bCi,*, at M = 2.0. 


NUMERICAL RESULTS AND DISCUSSIONS 


The vectors C,,/bC,,* and C,/bC,,* have been com- 
puted on the basis of Eqs. (23) and (31) for a range of 
reduced frequencies (0 < k, < 2.0) at two Mach 
Numbers (\/ = 10/7 and 2). Results of these compu 
tations are tabulated in Table 1 and are presented as 
the vector diagrams in Figs. 4 and 5, where the sub- 
scripts r and 7 have been used to denote the real and 
imaginary parts, respectively. 

From these vector diagrams, it is clearly seen that 
the imaginary component of the vector C,,/bC,,* is 
negative—i.e., the phase angle of the vector is nega- 
tive. On the other hand, for some low aspect ratio 
wings, the vector C,/bC,,* has positive imaginary com- 
ponent—i.e., its phase angle is positive. Furthermore, 
the magnitude of these vectors changes considerably 
with the aspect ratio. This is sufficient evidence that 
the wing-tip effects are of significant importance in 
purely rolling oscillations, especially for low aspect 
ratio wings. This finding is similar to the conclusions 
given in reference 1 about the vertical and pitching 
oscillations. 

The physical significance of the phase shift noticed 
above is most easily understood by critically examin- 
ing the possibility of instability in one degree of free- 
dom rolling oscillations of a rectangular wing. For 
this purpose, it is necessary to compute the work done 
by the oscillating wing system. For the rectangular 
wing shown in Fig. 1, the work done per cycle of roll- 
ing oscillation is 


W = Ff R, (dQ/dt), dt (45) 
0 


where FR is the rolling moment vector as defined in 
Eq. (2), dQ2/dt is the angular rolling velocity about the 
x-axis, and the subscript 7 is to denote the real compo- 
nent. From the given boundary condition, Eq. (1), it is 
clear that 

dQ OP 


= lim 


= —UAo exp (txt) (46) 
dt Y = 4+902 


By assuming a phase angle ¢ for the rolling moment 
vector, R may be written as 
R= R exp i(vt + €) (47) 


where |Rj represents the modulus of the R-vector. 
Thus, it is easily found that 


W = —(m/v)|R\AogU cos e = —N cose (48) 
where N is a positive quantity. Hence, for « < 
7/2, W < O-i.e., air absorbs energy from the wing 


system and, consequently, the system is aerodynami- 
cally stable. | > 7/2,W>O0 
supplies energy to the wing system and, then, there is 
Therefore, it may 


Only when (¢ i.€., air 
danger of aerodynamic instability. 
be concluded that, for purely rolling oscillations of a 
rectangular wing, there is possibility of instability when, 
and only when, the real component of the rolling mo 
ment vector becomes negative 

For the several cases being investigated [Figs. 4, 
5, and Eq. (44)], instability in roll would not be 
anticipated. 

A point of some interest which has also been borne 
out in the present note is the evaluation of the velocity 
potential in the mixed supersonic region by introducing 
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SCIENCES 


Values of (C,/bCz,*), and (C,/bC1,*): 


M = 10/7 


( C, ) ( +) 
ke bC1,* r bCz,” i 


——A.R.f = o — 
0 0.083333 0 
0.4 0.081729 — 0.007925 
0.8 0.077259 — 0.014473 
1.2 0.070854 —0.018601 
1.6 0.063790 — 0.019837 
2.0 0.057356 — 0.018357 
ke ——— AR. =3 —— 
0 0.047303 0 
0.4 0.046904 —0.001196 
0.8 0.045808 — 0.002023 
1.2 0.044268 —0.002201 
1.6 0.042636 —0.001615 
2.0 0.041262 —0.000318 
M = 2.00 
ke —— A.R.f = © ————_~ 
0 0.083333 0 
0.4 0.082513 — 0.004050 
0.8 0.080205 — 0.007427 
2 0.076833 — 0.009596 
1.6 0.073001 — 0.010262 
2.0 0.069365 — 0.009420 
k. ——- A.R. = 3 
0 0.060894 0 
0.4 0.060603 — 0.000477 
0.8 0.059810 —0.000661 
1.2 0.058700 — 0.000342 
1.6 0.057551 +0.000593 
2.0 0.056652 +0.002115 
ke , -A.R. = 
( 0.010416 
0.4 0.010707 
: 0.011647 
1.2 0.013134 
1.6 0.015120 


2.0 0.017526 


t Crp/bCr,* = (C,/bC1,*)A.R. = @. 





a correction potential term ®,; as in Eq. (8). This has 


(ici) 
bC1,*/; 
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(sci) 
bCi,*/ i 


———- A.R. = 5 ——_—— 
0.060508 0 
0.059688 — 0.003593 
0.057410 —0.006469 
0.054158 — 0.008084 
0.050602 —0.008194 
0.047414 —0.006892 
—— AR. = 1/8 
0.010416 0 
0.010668 0.003353 
0.011546 0.006555 
0.012870 0.009471 
0.014557 0.011995 
0.016473 0.014321 
——— AR. =5 

0.069471 0 
0.068983 —0.001817 
0.067627 —0.003203 
0.065675 —0.003823 
0.063524 —0.003500 
0.061601 —0.002248 
—— AR. =1 ——. 


0.027058 
0.027398 
0.028469 
0.030144 
0.032331 
0.034898 


1/p 
0 
0.004509 
0.009010 
0.013225 
0.017118 
0.020613 


0 

0.004165 
0.008206 
0.011904 
0.015159 
0.017889 
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Experimental and Theoretical Studies of 
Flame-Front Stability’ 


GEORGE H. MARKSTEIN? 


Cornell Aeronautical Laboratory, Inc. 


ABSTRACT 


Rich hydrocarbon flames burning in tubes were found to as- 
sume a cellular structure in the absence of turbulence in the ap- 
proach stream. The effects of type of fuel, mixture composition, 
and pressure on the phenomenon have been studied. A first- 
order perturbation treatment of flame-front stability gave results 
in qualitative agreement with experiments when the dependence 
of burning velocity on flame-front curvature was taken into ac- 


count. 
NOTATION 

a = acceleration 

ie; . Ae = constants 

rm = specific heat at constant pressure 
d = average cell size of cellular flames 
D; = diffusion coefficient of ith species 
F = aL/(S,°)? 

f(x), g(y, 2) = functions 

h = 2nr/X = wave number 

= V-1 

k = thermal conductivity 

3 = characteristic length of order of flame front 

thickness 

m = pu’S, = mass flow 

M = molecular weight of fuel 

M; = molecular weight of ith species 

ni = concentration of ith species 

Ni; = M;,(vi;/pDi) 

p = pressure 

Q; = heat released in jth reaction 

r = radial coordinate 

R = radius of curvature of flame front 
Su = burning velocity 

t = time 

T = absolute temperature 

u,v = velocity components 

W; = rate of jth reaction 

x,y = coordinates 

X = «/L 

Yo = T 

y; =ny,t=il,...,8# 

a = 6/uch 

B; = k/cypD; 

Y = a/h(S,°)? 

f) = stability parameter 

€ = pu? /p,° 
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angle between tangent of flame front and y- 


¢ = 
axis 

x = Lh = 2xL/xd 

Vij = number of molecules of ith species generated 
in the jth reaction (»;; is negative for those 
consumed) 

m = parameter that determines curvature depend- 
ence of burning velocity 

r = wave length 

Reina = wave length for maximum instability 

p = density 

o = uLlh 

T = parameter that determines curvature depend- 
ence of temperature of burned gas 

&n = coordinates of an element of flame front 

« = proportional to 

Subscripts 
u = unburned 
b = burned 
Superscripts 
0 = zero order 
1 = first order 
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INTRODUCTION 


—— STUDIES OF FLAME PROPAGATION in high- 
velocity gas streams! have shown that the inten- 
sity of turbulence in the combustion zone, computed 
from the observed flame speeds on the basis of existing 
theories of turbulent flame propagation,* * exceeded 
the level of turbulence in the approach stream by more 
than an order of magnitude. The high-intensity flow 
disturbances were assumed to originate in shear flow 
instabilities in the burned gas. 

However, based on experimental studies of the effects 
of artificial disturbances on laminar Bunsen burner 
flames, instability of the flame front itself as a source 
of flow disturbances has been suggested.‘ 


Subsequently, the observation of spontaneous break- 
down of rich propane flames into a number of separate 
cells® in the absence of turbulence in the approach 
stream provided direct evidence for instability of flame 
This phenomenon was later found also to oc- 
It appeared 


fronts. 
cur in many other hydrocarbon fuels. 
thus to be of sufficiently general significance to justify 
a detailed investigation. The experimental part of the 
present paper is primarily concerned with results of 
this study obtained with slow-burning mixtures of air, 
nitrogen, and hydrocarbons. 


Various authors have described phenomena, such as 
polyhedral*~ and truncated pyramid® burner flames, 
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Fic. 1. 


Setup for study of cellular flames. 


rising hemispherical flamelets,' filamentous and cor- 
rugated flame fronts,’ and distorted flame fronts in 
spherical bomb explosions,’ which seem closely related 
to cellular flame structure. The majority of these 
workers” ® ' and others'! favored an interpretation in 
terms of diffusion processes, while flow phenomena were 
considered to be only incidental. On the other hand, 
the only previous theoretical study of flame-front sta- 
bility’? was based exclusively on hydrodynamic reason- 
ing. In the theoretical part of the present paper, it is 
shown that, under certain simplifying assumptions, the 
physicochemical properties of the combustion zone 
may be introduced formally into the hydrodynamic 
treatment in the form of appropriate boundary condi- 
tions at the flame front. 


EXPERIMENTS 


Apparatus and Procedure 


The experiments were performed in vertical Pyrex 
glass tubes, open at the top and closed at the bottom 
except for the gas inlet. The flow velocity of the pre- 
mixed gases was chosen low enough to ensure laminar 
flow in the tube; the range of Reynolds Numbers was 
about 400 to 600. The mixture compositions were 
selected so that the flames remained almost stationary 





-MARCH, 1951 


within the tube, about 2 to 3 cm. below the top. In 
order to accomplish this over a wide range of fuel cop. 
centrations, the burning velocity was adjusted by add. 
ing suitable amounts of nitrogen. The proper composi- 
tions were found by starting with rich mixtures that 
gave diffusion flames burning on top of the tube and 
then gradually reducing the fuel concentration unti] 
the flames entered the tube. 

Average cell sizes of the cellular flames® were deter. 
mined by counting the number of cells over the whole 
tube cross section on flame photographs. Average cel] 
size was computed as tube inner diameter divided by the 
square root of the number of cells. The flame photo- 
graphs were taken under an oblique angle from below. 

In order to obtain good cell-size averages, it appeared 
desirable to take a fairly large number of photographs 
during each run, particularly since the flames were fre- 
quently distorted by downdrafts from the top of the 
tube. This was done by taking motion pictures instead 
of stills; good exposures were obtained with a rate of 
8 frames per sec. and f/2.3 on 35-mm. Eastman Lina- 
graph Pan film.* The experimental arrangement, 
consisting of a flame tube of 15.25-cm. inner diameter 
and the movie camera, is shown in Fig. 1. The shield 
on top of the tube served to protect the flame against 
drafts. The figure shows the setup in the Altitude 
Chamber for work at reduced pressures. Ignition at 
low pressures was accomplished by a pilot flame burn- 
ing from the nozzle visible on the right of the figure. 


Results 


With the exception of methane, all the hydrocarbons 
that were investigated gave cellular flames on the rich 
side of stoichiometric composition.f| The phenomenon 
was, in general, already visible in rich flames that 
burned on top of the tube as one or more approximately 
hemispherical cusps, convex toward the unburned gas, 
near the apex of the “inner cone.’’ As soon as the 
flames entered the tube, the whole primary combustion 
zone disintegrated into these cusps or cells. They were 
generally of fairly uniform size, convex toward the un- 
burned gas and separated by dark narrow ridges, and 
filled the tube cross section in a single horizontal layer. 
The cells were always in incessant irregular motion, 
which was too fast to be analyzed visually or on the 
movies taken at low frame rate. High-speed shadow- 
graph movies taken axially through the tube, which was, 
for this purpose, closed by a glass window at the bottom, 
revealed that the larger cells grew continuously at the 
expense of their smaller neighbors and split when they 
reached excessive size. 


Modifications of this structure were observed when 
the tube walls became heated in the course of consecu- 

* The author is indebted to F. Stratton for suggesting this 
photographic technique. 

+ With methane, extremely shallow cells were occasionally ob- 
served in lean mixtures. 
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STUDIES OF FL 


Moderately hot tube walls gave a structure 


tive runs. 
characterized by a ring of larger cells somewhere be- 
Extremely 


tween the wall and the center of the tube. 
hot walls gave a structure of a few large cells arranged 
in a symmetrical pattern and rotating around the tube 
axis. These phenomena have not yet been studied in 
detail. They were avoided during the determinations 
of average cell size by inserting cooling periods of at 
least 30 min. between runs. 

Cellular flame structures for various compositions 
and fuels differed not only in average cell size but also 
in what may be called their “‘strength.’’ This property 
can be described at present only qualitatively. In 
strong structures, the cells were almost hemispherical, 
completely separated by dark ridges, and in lively mo- 
tion. Weak structures had shallow cells, the separating 
ridges were often almost as luminous as the cells, and 
the motion was more sluggish. Cell size was more uni- 
form in weak structures. 

Extremely rich cellular flames, obtained without or 
with small amounts of added nitrogen, were character- 
ized by yellow radiation due to carbon formation, which 
was emitted partly by the outer mantle above the tube 
and partly in vertical streaks or bands that originated 
in the ridges between the cells. This phenomenon has 
also been observed by others’ and is analogous to the 
yellow radiation above the “‘open tip”’ of rich hydrocar- 
bon flames. Except with benzene (the only aromatic 
hydrocarbon studied), addition of moderate amounts of 
nitrogen and corresponding reduction of fuel concen- 
tration suppressed the yellow radiation completely. 
In benzene, however, it persisted almost until stoichio- 
metric composition was reached. 

The strength of the structures appeared to remain 
constant over a wide range of rich compositions. As 
stoichiometric compositions the 
strength decreased rather abruptly within a narrow 
range of fuel concentrations, until the structure dis- 
appeared completely and the flame assumed the shape 
The limit compositions at which cell 


were approached, 


of a smooth disc. 
structure disappeared were determined for all the fuels 
investigated except for the higher hydrocarbons that 
are liquid at normal pressure and temperature. It 
was found that the limit air-fuel ratios were independent 
of the amount of nitrogen added within the range over 
which the latter could be varied without causing flash- 
down or blowout of the flame. Limit fuel concentra- 
tions in terms of stoichiometric concentration are given 
in Table 1. 

The deviations from stoichiometric composition were 
within the errors of measurement. Thus, it seems 
plausible that the transitions took place exactly at 
stoichiometric composition. 

Contrary to a statement in the preliminary note,? it 
was found that, with each fuel, average cell size varied 
somewhat with composition. However, for fuel con- 
centrations within the range of about 1.20 to 1.40 times 
Stoichiometric, average cell size remained practically 
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TABLE 1 
Limit Compositions of Cellular Structure 


Fuel 
Concentration, 
Fraction of 


Hydrocarbon Stoichiometric 


Ethylene 0.983 
Ethane 1.030 
Propylene 0.970 
Propane 1.012 
Butadiene 1,3 1.019 
Butene 2 1.042 
Isobutylene 1.020 
N-butane 1.002 
Isobutane 0.997 


constant. This range, which also gave flame photo- 
graphs of best quality, was therefore chosen for com- 
paring various fuels and for studying the effect of pres- 
sure on average cell size. Beyond this range, cell size 
increased gradually toward richer compositions. In 
the case of n-butane, this increase was from 1.07 to 1.28 
cm. between 1.40 and 1.80 times stoichiometric fuel 
concentration. ‘ 

The appearance of cellular flames of various fuels 
in a tube of 10 cm. inside diameter at atmospheric 
These photographs were 
1.20 times 


pressure is shown in Fig. 2. 
taken with fuel concentrations of 
stoichiometric and about 30 per cent added nitrogen in 
They show that ethylene and ethane 
The flames of the 


about 


the atmosphere. 
gave comparably weak structures. 
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N- HEXANE 








1SO-OCTANE 


ISOBUTANE 


PROPANE 


Cellular flames in hydrocarbon-air-nitrogen mixtures at 
Tube inner diameter, 10.0 cm. 


Fic. 2. 
atmospheric pressure. 
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Fic. 3. Average cell size of hydrocarbon-air-nitrogen flames at atmospheric pressure vs. molecular weight of fuel. 
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Cellular flames in m-butane-air-nitrogen mixtures at 
Tube inner diameter, 15.25 cm. 


Fic. 4. (¢ 

various pressures. 
other fuels did not differ appreciably with respect to 
strength of structure. The photograph of the benzene 
flame clearly shows the streaks of yellow radiation, 
although it was taken with a dark-blue filter. 


Average cell size decreased markedly with increasing 
number of carbon atoms and decreased slightly with in- 
creasing number of hydrogen atoms in the fuel mole- 
cule: No difference between straight-chain and 
branched molecules was found. 
average cell size vs. molecular weight of fuel, Fig. 3,* 
showed that, with few exceptions, the data correlated 
remarkably well with the relation M'/d. = const. 
Excluding the values for ethylene, ethane, butadiene 
1,3, and benzene, M'“d values differed from their 
No theoretical explanation 


A logarithmic plot of 


mean by only 1.7 per cent. 
for this empirical result can be given at this time. 

The effect of reduced pressure on average cell size 
was studied by conducting the experiments in the Cor- 
nell Aeronautical Laboratory altitude chamber. JN- 
butane was selected as fuel. Fig. 4 shows the appear- 
ance of the flames in a tube of 15.25 cm. inside diame- 
ter at various pressures. A logarithmic plot of average 
cell size vs. pressure, Fig. 5, showed good agreement of 
the data obtained with tubes of 10.0 and 15.25 cm. in- 
side diameter and indicated good correlation with the 
relation p’“d = const. The deviations of the values 
for the lowest pressures were not surprising, since con- 
siderable difficulty to stabilize the flames and obtain 
good photographs was experienced at these pressures. 
Furthermore, because of the small numbers of cells at 
low pressures, wall effects may have falsified the re- 
sults. Significantly, the deviation started at a higher 
pressure with the 10.0-cm. tube than with the larger 
one. The significance of the result p’“d = const. will 
be discussed in the theoretical part. 

*In Figs. 3 and 5, the circles and crosses represent average 


values calculated from at least five photographs; the vertical 
lines indicate the spread of values from individual photographs 
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Fic. 5. Average cell size of n-butane-air-nitrogen flames vs. pressure. 


Exploratory experiments were performed with gas 
mixtures containing a fourth component in addition to 
hydrocarbon fuel, oxygen, and nitrogen. It was found 
that addition of small amounts of hydrogen to Jean 
mixtures of all hydrocarbons, including methane, caused 
strong cell formation. This is in agreement with Beh- 
rens’ findings.’ Ammonia caused formation of weak 
large cells in lean methane mixtures but caused an in- 
crease of cell size and decrease of strength if added to 
rich mixtures of higher hydrocarbons. Carbon mon- 
oxide had a similar effect on rich cellular flames. Use of 
carbon dioxide instead of nitrogen as diluent gave 
slightly smaller cells in butane-air mixtures. He- 
lum as diluent gas gave larger cells and _ sup- 
pressed cell formation completely if added in larger 
amounts. 

Studies of flame structure in fast-burning mixtures 
by means of high-speed shadowgraph and _ schlieren 
movie techniques are under way. A detailed discussion 
of their results is beyond the scope of this paper. The 
work done hitherto has indicated that there is a close 
connection between the phenomenon of vibratory 
flame motion! and cellular flame structure. Regard- 
less of the type of combustible mixture, cells were found 
to appear and disappear periodically in the rhythm of 


the oscillation during the vibratory stage. Fig. 6 is a 
sequence of high-speed shadowgraph movie frames, 
taken axially through the flame tube, which show this 
phenomenon. However, a significant effect of composi- 
tion was found: In those mixtures that were originally 
capable of forming cells, vibrations started at an earlier 
stage of flame travel and with a higher harmonic 
of the tube, while in originally 
tures, they started later and with the fundamental fre- 


noncellular mix- 


quency. 

It can hardly be doubted that the phenomenon of 
cell structure, observed to occur in the absence of 
stream turbulence, must be interpreted as an instability 
of the flame front. Flow phenomena undoubtedly 
play an important role; the resemblance to the phe- 
nomenon of cellular convection!" is striking. In anal- 
ogy to that phenomenon, one may regard each cell as a 
vortex ring, which, in the case of the flame front, how- 
ever, is superposed over the uniform flow through the 
flame front. In the following theoretical treatment, 
it is shown that at least qualitative agreement with ex- 
perimental results can be obtained by considering 
both the hydrodynamic and physicochemical aspects 


of the phenomenon. 








204 OF THE 


JOURNAL 





Fic. 6. 
motion in 10.0-cm. inside diameter tube. 
frames per sec. 


High-speed shadowgraph movie of vibratory flame 
Approximately 1,500 


THEORY 


Various problems of stability of phenomena that 
are governed by nonlinear partial differential equations 
have been attacked successfully by the method of 
first-order infinitesimal perturbations. Typical ex- 
amples of applications of this method are stability 
problems of parallel viscous flows'® and of cellular con- 
vection,'7—!9 

A rigorous formulation of the problem of stability of a 
plane flame front” must start with the Eulerian equa- 
tions, the equation of continuity, the equation of energy 
balance, and one diffusion equation for each chemical 
species participating in the reaction. Perturbations of 
velocity components, pressure, temperature, density, 
and concentrations of chemical species must be con- 
sidered. 

A difficulty peculiar to this problem arises from 
the 
calculation of the burning velocity of a steady plane 
been solved 


the fact that even the zero-order problem—i.e., 


has, in general, not yet 
This is primarily due to lack of knowl- 


flame front 
successfully. 
edge of the mechanisms of combustion reactions. 
Apart from mathematical difficulties, the presence of 
reaction-kinetic terms that form an integral part of the 
rigorous formulation thus constitute a serious obstacle 


for a successful solution. 
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Landau" evaded this difficulty by disregarding the 
physicochemical aspects of the problem completely. 
He regarded the flame front as a surface of discontiny. 
ity separating regions of constant steady-state velocity, 
density, and temperature. , 
and pressure only were considered in this purely hydro. 


Perturbations of velocity 
dynamic treatment. Thus, it is not surprising that 
the results obtained by Landau disagree with the ex. 
perimental results given above. While he found jn. 
stability for all flames (inasmuch as p, > pp» in all nor. 
mal flames), the experiments gave limit compositions 
Furthermore, he 
predicted a linear increase of instability with the wave 
number 27/A, whereas the observed uniformity of cell 


beyond which the flames were stable. 


size suggests a maximum of instability for a certain 
wave number. 

The treatment presented below occupies an inter. 
mediate position between a rigorous treatment and 
Landau’s. The relative simplicity of Landau’s ap- 
proach has been retained by restricting the validity of 
the treatment to wave lengths of the disturbance which 
are large compared to the thickness of the combustion 
zone. It appears permissible, then, to assume that the 
flow disturbance causes merely a bending of the reac- 
tion zone without affecting its internal structure di- 
rectly. 

Thus, the flame front may be still regarded as a sur- 
face of discontinuity, as far as the hydrodynamic part 
of the treatment is concerned. The physicochemical 
processes taking place within the reaction zone were, 
however, no longer disregarded; they were incorpor- 
ated into the hydrodynamic treatment by the boundary 
conditions at the flame front. A completely determined 
set of boundary conditions is obtained by specifying 
at each element of the flame front the local instantane- 
ous burning velocity and the temperature rise through 
the flame front. 
wave length of the disturbance, it was assumed that, at 


In view of the above restriction on the 


least in first order, local instantaneous burning velocity 
and temperature rise are functions of the local instan- 
taneous curvature of the flame front only. Thus, one 


may write 


ho 


I 


Si°[1 + w(L/R) + ...] (1) 
and 


T, = T,°(1 + 7(L/R) +... ] (2) 


where S, is the burning velocity and 7), is the tempera- 
ture of the burned gas, respectively, S,° and 7;° are 
their values for the plane flame front, R is the radius 
of curvature of the flame front, and L is a characteristic 
length of the order of the thickness of the combustion 
zone. For the purpose of the first-order perturbation 
treatment, only the linear terms in the expansions need 
be retained. 

The remaining difficult problem consists of calculat- 
ing the parameters w and r. Except under special as- 
sumptions, this cannot be done at present because of 
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lack of knowledge of the reaction mechanisms. In 
the Appendix, a general formulation of this problem is 
given, and its solution for a special case is indicated. 

F The hydrodynamic treatment, on the other hand, 
can now be carried out by introducing expressions (1) 
and (2) formally without further reference to chemical 
kinetics, diffusion, and heat conduction phenomena in 
the reaction zone. In order not to obscure the physical 
significance of the treatment by too lengthy deriva- 
tions, the further assumption 7 = 0 has been made in 
the following: This corresponds to the assumption of 
complete combustion, since, in this case, the tempera- 
ture of the burned gas must be independent of flame- 
front curvature. The more general case r ~ 0 offers 
no particular additional difficulties; its treatment will 
be published elsewhere. 

The treatment follows Landau except for the bound- 
ary conditions. Thus, the steady-state flame front 
is taken as the y-z plane, separating the region of un- 
burned gas x < O from the region of burned gas x > 0. 
The undisturbed flow velocities have the constant 
values u,° and u,° in the unburned and burned gases, 
respectively. Compressibility and viscosity are neg- 
lected. Thus, taking into account the continuity re- 
lation p,"u,° = py», the temperature ratio € across 
the flame front is equal to the velocity and inverse den- 
sity ratios 


= T;° 7 = ty” u,,° = py”, Pv° (3) 


Furthermore, in order to maintain a steady flame 
front, the conditions 
(4) 


u.° = $,°, m* = S,° = <,° 


must be fulfilled. 
across the flame front is given by 


p,? Du lds) (ec aia 1) 


Over this steady state, small time-dependent two- 
dimensional velocity perturbations u', v' and pressure 
perturbations p' are superposed. Neglecting second- 
order terms, the Eulerian and continuity equations for 
the perturbations become 


The steady-state pressure drop 


— pp = (5) 


Ou,,! 4,0 Ou,,! o 1 Op,’ 

ot "Ox p.° Ox 

SS ae ow (6u) 
Ot “On p.” OY 


(Ou,!/Ox) + (Ov,'/Ov) = 0 


for the unburned gas and an analogous system (6b) for 
the burned gas. 
usual procedure,'* '* a particular 


in y and exponential in ¢ is sought; 


Following the 
solution periodic 
thus, all perturbations depend on y and ¢ by 

g(y, t) = exp (thy + 6t) (7) 
where /; is the wave number, connected with the wave 
length \ by 
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Fic. 7. Kinematic conditions at an element of the flame front. 


h = 2n/d (7a) 


aud where 7 is the imaginary unit. The parameter 6 
plays the role of an eigenvalue of the problem. Sta- 
bility of the system is determined by the real part of 6. 
If any nontrivial solutions of Eqs. (6), subject to certain 
boundary conditions, correspond to positive values of 
the real part of 6, the system is unstable. 

Because Eqs. (6) have constant coefficients, the solu- 
tions are also exponential in x. By inserting solutions 
of the form of Eq. (7), it is easily seen that they must 
be of the form exp (Ax), exp (—hx), or exp [—(6/u°)x]. 

Of these solutions, only the first can be retained in the 
unburned gas and only the other two in the burned 
gas, since the perturbations must vanish at x = +o, 
Thus, with 


a = 6/u°h = 6/eS,°h (8) 


The solutions of Eq. (6) can be written 


Uy,'/Uy® = Arfi(x)g(y, t) ) 
V,'/iuy® = Arfi(x)g(y, t) 

Pu'/pu'(Uy®)? = — (ae + 1)Aifi(x)g(y, t) 
Uy'/Uy® = [Aofe(x) + Asfs(x) ]g(y, t) 
Up'/iuy® = [—Aofe(x) — aAafs(x) ]g(y, ¢) 
Po'/py?(tu®)? = (a — 1)Aofe(x)g(y, t) ) 


(9) 


where 

exp (x) 

exp (—hx) 
exp (—ahx) } 


filx) 
f(x) 
F(x) 


The A; are constants that have to be determined 
from the boundary conditions at the flame front. In 
order to establish the boundary conditions for the gen- 
eral case S, + S,°, one starts best with the rigorous 
kinematic conditions at an element of the flame front, 
Fig. 7. By definition of the burning velocity, the 
element at point £, 7 propagates relative to the unburned 
gas with a velocity S, normal to its surface and relative 


(9a) 


to the burned gas with velocity S, = Sy. 
Thus, the components of velocity of motion of the 
element £, 7 relative to the coordinate system are 


= uU, — S,cos¢g = &% — S,cos¢ 
v, + S, sing =%+ S sing 


dé/dt 


) 
(10 
dn/dt = f 
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where ¢ is the angle between the tangent to the flame 
front and the y-axis. 

In the first-order perturbation treatment, the dis- 
placements of the flame front are caused by the flow 
perturbations and are, thus, themselves small quanti- 
ties. Therefore, one may write, with sufficient accu- 
racy, 


(11) 


The displacement must depend on y and ¢ in the 
same way as the other perturbations. 


& = (1/h)Aag(y, t) 


cosg ~ 1, sing ~ g = O&/Ody 


(12) 


In Eq. (1) for the burning velocity, the radius of cur- 
vature may be approximated by* 


1/R = —07§/Ody? (13) 
The burning velocity is thus given by 
S, = S,°[1 + cAag(y, t)] (14) 
with 
o = wLlh (15) 


Finally, the pressure drop across the flame front is 
given by 


Pu — Po = puS.2(e — 1) (16) 


Inserting Eq. (9) for x = 0 and Egs. (11), (12), and 
(14) into the boundary conditions (10) and (16), one 
obtains, after canceling zero-order terms on the basis 
of Eqs. (4) and (5) and neglecting second-order terms, 
the following system of linear homogeneous equations 
for the constants A, As, A3, and Aa: 


A; — (ae + o)Ay = 0 
A, + A; — ea + o)As = 0 ( 
A, + Ao + aA; — (€ — 1)A, = 0 

(ae + 1)A; + (a — 1)A2 + 2o(e — 1)As = 0 | 


The system has a nontrivial solution only if the con- 
dition 


) 


(17) 


‘1 0 0 —(ae + o)!| 
0 1 1 -ea+a)| _ ' 
1 L s—aae* 
la +l a—1 0 = 2o(e — 1)| 


is fulfilled. This is a cubic equation for a with the roots 


1 
=e" 5 +: 
l Pa 
| -« +o) += (l+te—-+e0? — 200)"*|h (19) 
€ 
a3 = 1 


The root a3 must be discarded, since it leads to van- 
ishing solutions (9). The negative sign of the square 
root in a, 2 corresponds to a stable solution. 


* Thus, curvature is defined as positive when the flame front 
is convex toward the burned gas. 
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Thus, stability will be determined by the solutioy 


S,° 
; xe x 
Lite 
P 1 1/9 
{[1+.-! + (ux)? — 2anx| ~ (1 +ux)f (20) 


where 


2rL, A (21) 


and Eqs. (8) and (15) have been taken into account, 
This solution is plotted in Fig. 8 for « = 5. Lan. 
dau’s solution corresponds to » = 0. It is seen that 
the effect of curvature on burning velocity is stabilizing 
if » > O—1.e., if the burning velocity is larger when the 
flame front is concave toward the unburned gas. This 
is the type of dependence one would normally expect 
on the basis of simple thermal theories of flame propa- 
gation. The stabilizing effect increases with increas- 
ing wave number, and, thus, a maximum of instability 
for a certain wave number is obtained in qualitative 
agreement with the experimental results. 
in the Appendix that for the special case 


D,; = k/cyp 


It is shown 


C= 1...) (22) 


where D, is the diffusion coefficient of the 7th chemical 
species, k is the thermal conductivity, and Cc, is the 
specific heat at constant pressure; one obtains u = | 
with 
L = k/cpp.°S,° (23) 

Fig. 8 shows that, for this case, maximum instability 
occurs for 27L/X,,q,, =. 0.2. With the reasonable 
value L = 0.03 cm. at atmospheric pressure, a wave 
length of about 1 cm. is obtained in good order-of- 
magnitude agreement with the experimental average 
cell sizes. 

As to the dependence of average cell size on pressure, 
under the assumption that e, u, and k/c, are independ- 
ent of pressure, the theory predicts 


Amar. & L & (py°Sy°)-! & (pS,°)— 


maz- 
The experimental result d « p~‘* would thus yield 
the relation S,° « p~'/'. The same relation has been 
theoretically predicted*! and experimentally verified 
for various hydrocarbons.”? However, work on the 
pressure dependence of the burning velocity of butane- 
air mixtures** did not substantiate this relation and 
gave results that depended critically on fuel concentra- 
tion. Further experimental work on burning velocities 
at reduced pressure is thus needed before the result d « 
p~‘*’* can be interpreted without ambiguity. 

Eq. (20) still does not explain the existence of com- 
pletely stable flames, since a range of instability exists 
for any value of u. However, for large positive values 


+ This maximum is actually a saddle point for the real part of 
5 in the complex x plane. A discussion of the physical signifi- 
cance of this resuit will be published elsewhere. 
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of u, this range could easily correspond to cell sizes 
larger than the tube diameters used in the experiments, 
and instability could thus have escaped observation. 
Furthermore, the effect of gravity has been neglected 
sofarin the treatment. The influence of gravity acting 
in the direction from the burned to the unburned gas is 
equivalent to an acceleration a in the opposite direc- 
tion. It causes an additional pressure drop across the 
flame front proportional to the displacement £ and the 
difference of densities. Thus, 


(p,! ye Po') 
This gives an additional term in the last of Eqs. 
(17), which becomes 


(ae + 1)A; + (a — 1)A2 + (€ — 1) X 
[20 + (y/e)|As = 0 


= a(px" — pp )E (24) 


acceleration 


(17a) 
where 
a/h(S,°)? 


y= (25) 


The final result is 


S.° l 
Afrten i+ 


(: — *) *)" —(1+ wo} (20a) 


xe 


2 «a 9 ~_ 
Tiae UX) 2eux 


Plot of stability parameter vs. 


wave number for several values of u 


where 


F = aL/(S,°)? (25a) 


Fig. 9 shows a plot of this expression for e = 5 and 
uw = 1. It shows that the effect of gravity is stabilizing 
and increases with the wave length. Thus, a finite 
range of wave lengths for instability is now obtained, 
and, for sufficiently large values of yu, stability over the 
whole range is possible. It is interesting to note that the 
acceleration term in Eq. (20a) is identical with the ex- 
pression for the effect of acceleration on the stability 
of an interface between two liquids.” 

Finally, under conditions of vibratory flame motion, 
alternating accelerations are superposed over the effect 
of gravity. Thus, instability would increase and de- 
crease periodically corresponding to the family of 
curves, Fig. 9. The results of the high-speed movie 
studies, Fig. 6, are in agreement with this conclusion. 


APPENDIX 


It is desired to calculate the dependence of burning 
velocity and temperature of the burned gas on curva- 
ture of the flame front. By restricting the treatment 
to small values of L/R, it can be formulated as an eigen- 
value perturbation problem. 
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The basic equations for a plane steady flame front are Introducing the parameters 
d’T° ar wu? = (cym°L)/k, Ni = Qi/k 
5 i Cpm® + Vv k (27) 
dx° dx Ny = M,;—, B= emul...) "7 
pD,; CppD, 
Ol T?, HY; 6.65%) = 0 
F | and replacing the independent variable by 
d*n,° dn;,® t (26) A = 2/L (28) 
pD; — — 0 _f. 
dx? dx | ; ie 1s by 
and the dependent ones by 
M;, v,jw;(7, n,°, «eg n,.*) = 0 Y,° = 7 Y° =n, («= 1 n) (29) 
. , i t ’ ee , aaa 
j . 
me? = 9.8 J Eqs. (26) can be written 
ies . . ; ay, aY;? : : : 
where ,° is the concentration of the 7th chemical spe- — — Buy® —- + N,,w;(Yo°, ..., Yn°) =0 
cies, Q; is the heat released in jth reaction, w,is the rate @*” dX 
of the jth reaction, »;; is the number of molecules of (¢ =0,...,m) (26a) 


ith species generated in the jth reaction, and ;, is the 
molecular weight of the ith species. The heat con- 
ductivity k and the product pD; have been assumed to 
be constant.* The mass flow m®, which determines 
the burning velocity, is an eigenvalue of the nonlinear 
system (26). 

* Eqs. (26) lack generality in other respects.24 They are used 
in this form here rather to show the method of attack than to 
lead to a rigorous result. 


The general equations for a cylindrically curved 
flame front are 


(= l =) dT 
k — + — cm + 
dr? 


r dr dr 
) QO)w;(T, m, ..., Wn) = 0 
j 


(and similar equations for the n;). 


(30) 
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If the range of the radial coordinate r over which 
Eqs. (30) have to be integrated is small compared to 
average value R, it is permissible to replace 1/r by 


its 
Furthermore, the dependent 


|/Rin the second term. 
yariables and the eigenvalue may be developed in terms 
of L/R and the variable r replaced by X = (r — R)/L. 

If this is carried through, canceling of zero-order 
terms and neglecting higher-order terms in L/R leads 


to the system 


ry; AV ‘ 
xe" dX » il 
+? bey Vi = (aut — 1)" (an) 
‘= AT le — <3 
mt \OVi) y= vs 7 dX 


where the reaction rates have been developed into 
Taylor series with respect to the zero-order values and 
where substitutions analogous to Eqs. (27) and (29) 
have been performed. 

Assuming that the solution of Eqs. (26) is known, the 
problem consists of finding a solution to Eqs. (31) 
which fulfills the proper boundary conditions. jy! is 
the first-order perturbation of the eigenvalue and is 
thus the parameter that determines the dependence of 
burning velocity on curvature. The value that Y%! 
assumes at the upper boundary determines the depend- 
ence of temperature of the burned gas on curvature. 

Solution of this system is, in general, a difficult prob- 
lem. For the special case of complete combustion, a 
solution can be found if, in addition, the conditions 


Bom 1 (¢ = 1, ...,%) (32) 


In this case, the right sides of Eqs. (31) 
vanish for u' = 1. By differentiating Eqs. (26a) with 
respect to XY, it is seen that the derivatives 


are fulfilled. 


Zi° = dY;°/dX (33) 
fulfill the equations 


bi dZ,° 4 N (5) a 
45 2D . 23 OV ynve 
J 1 


(34) 


d 27 0 
dX? 


which are formally identical with the left sides of 
Rgs. (31). The Z;° vanish at the boundaries, since the 
Y,° are constant in the unburned and burned gases. 
In the case of complete combustion, temperature and 
concentrations in the burned gas must be independent 


of curvature; thus, the Y;! also vanish at the bound- 


aries, 
Thus Y;! « Z;° is a valid solution of Eqs. (31) in the 
special case 8; = 1 with the eigenvalue uv! = 1. By 


selecting the value of Z so that »° = 1, this solution 


becomes 
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L k 
De = See l see L —= (35 
| + R + | ( ==) 35) 


(32a) 


for 


D, = k/c,p («= 1,...,m) 


This result is physically reasonable, since under con- 
ditions (32a) the concentration and temperature dis- 
tributions will not be displaced relative to each other 
by the effects of curvature on heat conduction and dif- 
fusion. Thus, under these conditions the effect of cur- 
vature is independent of the reaction kinetics. 
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Thickness of a Steady Shock Wave 


J. J. Bernard , 
Office National d’Etudes et de Recherches Aéronautiques, Paris 
November 10, 1950 


Us THE ABOVE TITLE (JOURNAL OF THE AERONAUTICAL 
* ScreNcEs, Vol. 17, No. 9) R. von Mises gives some results 

extending the solutions of M. Morduchow and P. A. Libby! who 
considered the case of Prandtl Number 3/4. This case had al- 
ready been studied in France some years ago by M. Roy,? and I 
should like to point out some results obtained recently* by gen- 
eralizing Roy’s method. 

The symbols are those defined in reference 1: V 
& = x/k\y; @ = T/T); M = Mach Number at minus infinity; 
P = Prandtl Number; and » = pof(0) gives the law of variation of 
the viscosity. 

Thus, the basic equations are: 


/ 
= u/uo; 


4 iV 
0+ 7OV' — (1 + oY - = YMfO)V = “0 @ 
0 1 P 7 ¥ f(0) de 
=- = yay? 1 OS a ae a 
+“) 9° le ale s~-i MP a 
Y =a 
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(1) It is easily seen that the law f(@) does not affect the ex. 
pression of @ as a function of V, since this function 6(V) is a solu. 
tion of the following functional equation, derived from Eqs. (1 
and (2): 


4P\ . . 4P y¥-1 
é+{f1—- 3 & 6e-*dz = | (ie —_ M2 (1 — v9] 
« o ~~ 


3) | 
where . 
. aoa | 0G 
3 = yM*S V dV/[0 + yM? + (1 + yM2)V | 
as f(0) is eliminated, Eq. (3) can be integrated by an iteration 
method. The complete solution of the problem is then possible | 
(2) In the particular case of f(@) 1, the gradient (dV dt) is 
given, asa function y(V), by the equation 
: ki As as 
Vyy’ = M?P(1 — V) (V — a) + 
8 ¥ 
P 3h... 3 y¥+1 
uf ( + ) —-(1+ a) | —y > (4 
Y 2 8 ¥ : 
where 
t 
a= [24+(7- 1) M?] (vy + 1) M? 
The complete solution yy) may be obtained for P = 3/4; for any 
other value of ?, an approximate solution is obtained by writing 
y = y + hand linearizing; thus, 
3 y¥ —1(1 — V)&V — a) 
h(V) = M(=-—P “x aq 
4 ¥ J 
7] | 
/ (1 —a) —(V at} F I ( V — adt + alta 
0 
where 


B = [4P + 3(7 — 1))M?/3(M? — 1) 


By means of Eq. (5) the gradient (dV/dt) = y(V) is determined, 
at least for P Y 3/4; this solution seems acceptable even for large 


values of P. Eq. (1) gives then @(V) for f(@) = 1. This expres- 





° . ° ‘ pe 
sion of 6 remains unchanged for any function /(@) and thus is ob- ‘YI 
é ? : ee coup 
tained more easily than by integration of Eq. (3). h 
i ‘ e ‘ ° as 
(3) Shock-wave thickness / against Prandtl Number is plotted 
. # ; ee ‘ Pe deco! 
in Fig. 1 for f@ = 6. With Roy’s definition, shock-wave tt 
, : : ; ne s 
thickness is the distance between the two points where the great- oe 
: cnn : i or t 
est curvature of velocity distribution occurs (definition of second : ; 
t 
order). 
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FIG. 1. 


type where the collective blade pitch angle is kinematically 
coupled the coning angle. This “pitch-cone change” 
has the following effect: (a) The slope of the a; — @ curves 
becomes negative—i.e. the rotor is statically stable (b) 
the slope of the cr — a@ curves is appreciably smaller than that 
for the ordinary fixed-pitch rotor. 

It may be of interest to mention that these aerodynamic 
characteristics are similar to those of an ordinary fixed-pitch 
The curves show the thrust 


with 


and 


rotor in autorotation, see Fig. 1. 
coefficient cr (above) and the tip path plane inclination a, 
(below) against the attitude of control plane a. The full lines 
refer to a rotor with pitch-cone change and are taken from the 
paper by Hohenemser. The broken lines correspond to wind- 
tunnel tests with an ordinary autorotating rotor with 4° pitch 
setting and 0.06 solidity, where the individual values marked 
along the broken lines indicate the advance ratio uy. 

Hohenemser has shown that a helicopter with pitch-cone 
change is dynamically stable if the moment of inertia of the 
fuselage is sufficiently small or, more accurately, if the non- 
dimensional quantity J/m is smaller than the boundary value 
given in Fig. 13 of the original paper. Apart from this stability 
bovndary, however, nothing has been said about the grade of 
Stability or instability which may be expected. 

To answer this question, some numerical investigations have 
been carried out. The calculations are based on the assump- 
tions of Hohenemser’s paper and refer to some well-known 
helicopter types such as Bell 47D, Sikorsky S-51, and Sikorsky 
R-4B. It has been assumed that the original rotor of these 
helicopters has been replaced by a rotor with pitch-cone change; 
the effect of the fuselage on the dynamic stability has been 
neglected. The investigations show that, in these examples, 
the quantity 7/m lies in the range of approximately 0.06—0.08. 
According to Fig. 13 of the subject paper, this means that, in the 
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condition, the types mentioned above are 
In the case 


“constant torque’ 
dynamically unstable throughout the speed range. 
“constant rotor speed,”’ the most favorable type (Bell 47D) 
is dynamically unstable up to 0.3 advance ratio and is dynami- 
cally stable above this figure. For the Sikorsky R-4B at 0.35 
advance ratio, the damping coefficient of the unstable oscillation 
was calculated. It was found that, in the ‘‘constant-torque”’ 
condition, the time to double the amplitude of a disturbance 
increases from approximately 5 sec. for the ordinary rotor to 7 
sec. for the rotor with pitch-cone change. 

These figures indicate that the improvement 
longitudinal stability (stick fixed) due to pitch-cone change is 
probably of minor importance. This applies to both the single 
rotor and the laterally arranged two-rotor configuration. Of 
course, it does not exclude the possibility that the general 
handling characteristics of a helicopter with pitch-cone change 
show a more appreciable improvement. 


in dynamic 


+ 
Buckling of a Pretwisted Pin-Ended Column 


G. V. R. Rao 

Department of Aeronautical Engineering, New York University, 
New York, N.Y 

December 1, 1950 


or STRUCTURAL MEMBERS under compressive loads are 
supported at the ends by pins. It also often happens that 
the pins, because of misalignment, may have to be twisted to 
Thus, the strut is subjected 
The stability of 


make them parallel to one another. 
to a torque in addition to the compressive force. 
such a strut is investigated, and it is found that the effect of the 
torque is to reduce the buckling load below the Euler load. 
Gregory! has considered the same problem, but he used incorrect 
end conditions and obtained a buckling load higher than the 
Euler load of the column. 

The coordinate system and the end forces are shown in Fig. 1. 
With frictionless pins, the force R in the direction of the pins 
must be taken as zero, and the force S in the perpendicular 
direction evidently is not zero. The differential equations will 
then be 
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— EI (d*x/ds?) = Px + T(dy/dz) — M 


Py — T(dx/dz) — Sz } 


— EI (d*y/dz*) = 
The following general solutions: 
x = A,cosnz — Asin 1s + A; cos mz — Aqysin rez + 
(M/P) — (ST/P?) 
y = A, sin mz + Aocos nz + Az sin res + Ay cos rez + (S/P)zs 


will satisfy the differential equations [Eqs. (1)], provided 


1 0 1 
0 l QO 
0 = 0 
cos p; —sin Bp, cos Be 
sin p; cos p; sin By 
—B, sin B; — Bi cos B; — Bz sin Be 


where 8; = /r,; and B2 = Ire. 


The determinant can be simplified to yield the following 
conditions: 
Bz sin Bo(1 — cos B,) = B; sin B,(1 — cos Be) (5) 
Thus, either 8; = Bs, which leads to a trivial solution, or 6; = 7 
and 6 = — mw—i.e., % — tr = 2r/l. Substitution of this in 
the condition (2) yields 
(T/EI)? + (4P/EI) = 4n?/I? (6) 


This equation is the same as that given by Timoshenko? for the 
buckling of a column subjected to a twist and thrust with uni- 
versal joints at its ends. The result might seem incongruous, 
but, if we examine the bent shape of the column in half its 


length, the equality of the two cases becomes apparent. When 


On the Equations of Longitudinal Stability—I! 


John W. Miles 
Associate Professor of Engineering, University of California at Los 
Angeles 


December 22, 1950 


D* NEUMARK (R.A.E., FARNBOROUGH, ENGLAND) has 
pointed out that the result E = 0 in Eq. (14e) of a recent 
note! on this subject is at variance with the usual result for free, 
power-off flight. 
X0) = Xup = Xq, Whereas we can state correctly only xg, = xq, and 
Eqs. (14) should read 


This was due to an incorrect assumption that 


A=] (14a) 
B = (Sq + mg + Me) (14b) 
, zm 

C = uma + a ”) (14c) 


xm som zm zm 
ee igi 2 ails CE | es 4 he | 
uo a Uy a g Uy 8, ui a 
(14d) 
‘ . ' Zz m\ 
E = p (xa — Xu) ae (14e) 


the notation being that of reference 1. 

Further, the derivatives with respect to n,, representing virtual 
inertia in the direction of flight, may be established as small by 
direct dimensional considerations, and m,, generally is small, 
vanishing for trimmed flight at low speeds. With these approxi- 
mations, Eq. (14d) reduces to 
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—EIr?+Tr+P=0 " 


The unknown constants A;, A», A;, As, M/P, and S/P can} 
solved for by substituting the boundary conditions, 


x = y = dx/dz = 0 9 
atz =Qand/. Asaresult, x = y = 0 everywhere, or the follow. 
ng determinant is zero: 

0 1 —T/P 
1 0 0 
= 0 0 
- . 0 j 
—sin Bo ] —T/P 
cos By 0 l 
— Bz cos Be 0 0 


the ends are universally jointed, the bent shape of the columns 
a complete reflection about the centerpoint, and in the present 
case it is a reflection in a central plane. 

Eq. (6) can be modified as 


P = (nr*EI/l?) — (T?/4EI) 
which shows that the buckling load is below the Euler load of the 


column. 
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Stability, p. 167; McGraw-Hill Book 


Company, Inc., 


D= 


x m xsm 
7 + + Su (X09 — Xwo) (mg + me) 
uy @ uo aq 


(14d 


Hence, the conclusion that m4 is the only unsteady flow term that 
need be included in the longitudinal stability equations remains 
unchanged 
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Note on the Relation of Lifting-Line Theory to 
Lifting-Surface Theory 


Eric Reissner 
Department of Mathematics, Massachusetts Institute of Technology 
December 20, 1950 


oo FOLLOWING IS CONCERNED WITH linearized lifting-surface 
theory for the pressure distribution on wings in subsonic 
motion and its relation to lifting-line theory for the spanwise 
variation of lift. It is shown that in certain cases lifting-line 
theory may be obtained by means of systematic developments 
with regard to aspect ratio from lifting-surface theory. It is also 
indicated by means of a typical example that such developments 
are not possible in all cases where lifting-line theory is considered 
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_* TABLE | 
= s 0 0.1 0.2 0.5 1.0 2.0 3.0 4.0 5.0 
0 0.196 0.32 0.58 0.84 1.11 1.23 1.31 1.35 
Fe 0 0.196 0.32 0.58 0.81 0.92 


to furnish an approximation to the results of lifting-surface the- 


ory. 
Consider a lifting-surface of infinite span, of uniform chord 
9b, moving with subsonic velocity U in a fluid of undisturbed 
density p. The integral equation of the linearized theory for this 
problem may be written in the form 


a 1 


1 a 1 . 
al + Ke — by) ] dean (1) 
2r ’ -; OnLy— 1 


K= V(x — §&£)? + (y — 0)? — ly — 0] (2) 
(x — &)(y — ») 
The dimensionless coordinates x and y are related to the physical 


coordinates XY and Y in the form 

x = X/b, y = (Y/b)V 1 — Mt? (3) 
and the dimensionless downwash function w is expressed in terms 
of the local lifting-surface slope a( .Y, }) by 


w = (2pU2/V1 — M?)a (4) 


The function p gives the pressure jump across the lifting surface, 
M is the Mach Number corresponding to the velocity U, and 
principal values are to be taken of the integrals in Eq. (1). 
Assume a downwash function w of the form 
w(x, y) = W(x) cos (y/s) (5) 
The aspect ratio of a portion of the lifting surface between con- 
secutive nodal lines of the downwash distribution follows from 


Eqs. (5) and (1) as 
A.R. = rs/2V 1 — M? (6) 


The double integral equation [Eq. (1)] may be reduced to a single 


integral equation by setting 


P(x, y) = 
Combining Eqs. (7) and (1) requires evaluation of two infinite 
integrals with respect to 7. With a new variable of integration 
st, these integrals may be written as 


So ee y : 
sin ¢ = — 7 Ccos (8) 
is s 


a . 
sin- A(x — £, 9 — n) dn = 
§ 


y ear St ae : 
— cos sin ¢ K tr & (9) 
Ss Ss 


Setting as an abbreviation 


P(x) cos (y/s) (7) 


{, defined by y — 7 = 
follows: 
’ _  «< dn y 
sin = 

s 77> © Ss 


— sin ¢ K( 


pa 0 


F(s) = z|, ¢) dg (10) 


the resultant one-dimensional integral equation becomes 
1 P(é) 


1 1 , 2 x—¢ 
P(é)| 1+ - F dt 
7 x—é& 2s J_s T s 


d§ + - 
—] 
(11) 


The function F may be expressed in terms of tabulated func- 


W(x) = 


tions. From 


: Vite—ty 


sin ¢ $ 
0 aa § 


F(z 


: Vili+uw—u 


= sin |clu du 
- 0 u 


(12) 


follows by appropriate differentiation and integration by parts, 
fors > 0, 


d q dF i cos fu ’ 
F+sz = ; du = Koc 
dz dz So V1 + u? 


In Eq. (183) Ao isa modified Bessel function of the second kind. 
From this there is obtained, by integration and by making use of 
some properties of the modified Bessel functions, the representa- 


40% 


For sufficiently small values of c, this implies the approxima- 


(13) 


tion 


l 
F(z) = Ko(s)ds + Ki(s.) — | (14) 


tion 


F(z) ~ [0.808 — (1/2)In |z| | s (15) 


‘.(s) = 


which is based on the leading terms in the power series expan- 
sions for Ky and Aj. 
For large values of s one has the limiting relation 


(16) 


F( +o) = + 4/2 


Table 1 contains numerical values of F and F,. 
It may be noted that the approximation (15) is certainly ade- 
quate as long as ¢ is less than 1, and this means that as long as 


jx — £|/s is less than 1; or, since |x — &| < 2, as long as A.R. > 
x/V1 — M?, the integral equation [Eq. (11)] may be written 


in the form 

71 

> P(&) dé + 

x—-—¢ 2s 1 
1 


1 i Ins In|jx -—¢€ 
P(é) (x — £){ 0.808 + — = : dt 


ws? = 


(17) 


Eq. (17) is a special case of a class of integral equations which 
Eichler! has shown to be explicitly solvable 
The following observations may be made: Omission of the 
last of the three integrals on the right of Eq. (17) leaves an equa- 
tion the solution of which corresponds to Prandtl’s lifting-line 





theory 
A development of the solution of the complete equation [Eq. 
(17)| for sufficiently large values of s will be of the form 
l In s 1 
P= PP, +-P, + — 2. +— Ps + (18) 
Ss a $$” 


The term Po is the exact solution of the two-dimensional theory. 
The term P, is correctly taken account of by lifting-line theory. 
The term P, is not taken account of by lifting-line theory, since 
this theory does not introduce logarithmic terms in s. 

In order to obtain a result corresponding to Weissinger’s 
improvement of lifting-line theory * we may proceed as follows: 
Write 

“1 
(19) 
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TABLE 2 = 
5 0 0.2 0.25 0.33 0.5 1.0 2.0 5.0 10.0 x 7 
L/lp> | 0.5 0.565— 0.58 0.605 0.65 0.755 0.86 0.955 0.985 1.0 = 





and integrate Eq. (11) in the form 


1 —— 1 1 m 
; 1+x 1 [1 +x dx 

W (x), /—— dxk& = ~ P(é) - 
ii ” Vi — <£ ; cs - cof NG —“2L£— 


x—e 

fret feel PE] 
P(é) « i ere > dé 

taf, (€) ~1VWil—=« s f 
(20) 

In Eq. (20) 
1 | “ 1 | 
1+x dé ji + x 

= 7, dx = 21 

¢ Vs, " Lv * 


and in the term multiplied with F the function P is approximated 
by 

P(E) = (L/w) V (1 — £)/(1 + £) (22) 
There remains a double integral over F which is approximated as 
follows :3 


{ff i+ x po *) cag = = #(*) 
= «a Wi-sWVite s s 
(23) 


Eq. (20) now reads 


: 
l+x Boas 

W(x) | dx = z{1 er sehe r( } 24 

J. (x) V1 = *# a” TS s a) 


In view of Eq. (15), we have, for sufficiently large values of s, 


: 1 1 
W(x) jl + 2 de ro 7, E + = + (0.808 + 3 + - , In ‘| 
—1 V 1—<x 2s s? 
(25) 


and thus, in this improvement of lifting-line theory, we do have 
a logarithmic term in s, as in the solution of the lifting-surface 
theory problem. However, it may be observed that the numeri- 
cal coefficient of the In s-term in Eq. (25) will not, in general, 
agree with the corresponding coefficient following from lifting- 
surface theory, since, in lifting-surface theory as formulated in 
Eq. (17), the value of Z depends on the function W(x) in a less 
simple manner than predicted by Eq. (24). In this regard, an 
improved approximation will be obtained by using as in reference 
3, instead of one equation for L, two simultaneous equations for L 


and 
1 
w= f P(€) & dé 
—1 


The result (25) for sufficiently large values of s may be con- 
trasted to the result for small values of s. In view of Eq. (16), we 


have, in this case, 

& ; T T i] wv 

We) Ved ~w Lili t+—+—$ 4L- 

x 2s 2s f Ss 

and, consequently, the value of Z in this case is one-half of what 

it is according to the original lifting-line theory. This result is in 
accordance with a recent observation by Lawrence.‘ 

If we denote the value of LZ according to lifting-line theory by 

Lp, we have, for the ratio of L of Eq ca to Ly, 


Lb 1 + (x/2 


(26) 





Lp 1+ (x/2s) + (1) s F(1/s) 

Table 2 contains numerical values of the ratio L/L» as function 
of s and, in view of Eq. (6), as function of A.R. and M. 

While it is seen to be possible to obtain lifting-line theory for 

the infinite lifting surface with spanwise periodic downwash dis- 

tribution by means of systematic developments in terms of the 


(27) 


aspect ratio parameter s, it may also be seen that ac orresponding 
procedure is not possible for a finite lifting-surface with finite tip 
chord. We need only consider the semi-infinite lifting surface o 
constant chord, with constant downwash function w(x, ¥) for 
y 2 0. In this case the aspect ratio is infinite. Lifting-ling 
theory may be used to obtain an approximate result for the lift 
distribution that is zero at the tip y = 0 and approaches a finite 
limiting value for increasing values of y. This distribution wil}, 
however, be different from the distribution obtained by solving 


the double integral equation [Eq. (1)] with the integral J_°, re. 
placed by the integral J, : 
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Theodorsen! 
dimensional stability characteristics by an 
transformation and finite aspect ratio correction. 

In airplane dynamics, a system of axes moving with the body is 


THEORY of 
three- 


NONSTATIONARY AIRFOIL 
the determination of 


TWO-DIMENSIONAL 
is applied to 
appropriate axis 


commonly used, while in reference 1 Theodorsen uses axes parallel 
and perpendicular to the average motion of the wing. From 





Fig. 1 it is seen that the relation between velocities in these 
coordinate systems is 
h = weosa — usina 
v = ucosa+wsina 
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where h, v, and a@, as indicated in Fig. 1, are those used in ref- 
erence 1, u and w are velocities along the stability axes, and 9 is 
the angle made by the airfoil axis with some fixed reference. 


If ais a small angle and w < wand h < v, then, approximately, 
— ua, v= t& 


To the same approximation P = Z, 
perpendicular to the average wind and Z is the normal force 
Also, the moment remains un- 


h=w 


where P is the lift force 


referred to the stability axes. 
altered by the transformation, so it is only necessary to express 
Zand M in terms of the new velocity coordinates w and @. 
From Fig. 1, it is seen that 6 = a + c, where c is a constant that 
may be made zero by taking the reference axis parallel to the 
average velocity. Also, using A for the ordinary instantaneous 
angle of attack, 


(h u) + 6 


From reference 1 (Eqs. XVIII and XX, omitting the 8 terms, 


A= 


w/t = 


Z=P = —nxpb?(va + h — baa) — 
2xpvbC(k) {va + h + b[(1/2) — ala} (1) 
M= -xpb*| (3 _ a)obe + (2 + aa _ abi |+ 
3 ; 1 , 
2xpob*( a + \ (| va +ht+ (1 - a)é | (2) 


Substituting h = w — uO and a = 6 gives 


Z = Z(w) + Z(8) (3) 
M = M(w) + M(@) (4) 


where the separate parts due to the independent coordinates w 


and @, in analogy with the stationary stability derivatives, are 


Z(w) = —xpb%w — 2irpubC(k)w (5) 

Z(6) = xpb3a6 — 2rpub*C(k)[(1/2) — aj6 (6) 

M(w) = xpb2abw + 2mpb2ula + (1/2)]C(k)w (7) 
1 . 1 a » 

M(6) = _se{ (2 _ a) ubé a »(2 oe a?) + abué | + 


2mpub*(a + 1Vccwo(2 ~~ a)é (2 


In the equations above, C(k) depends on the distribution of 
vorticity in the wake and therefore on the history of airfoil motion 
Thus, while it appears that w, w, 6, and 6 are independent, 
C(k) can only be evaluated for a particular history for each of w 
and @(i.e.,hand a). Selecting wand @ for independent variables, 
all the time derivatives w, w, ..., and @, 6, 
through the parameter time; therefore the force and moment 
are determined by the complete motion. They can be separated 


., are determined 


4 


FORUM 215 


into parts due to each derivative of the motion only when it is 
slow enough for the circulation to be approximated by its quasi- 
steady value. 

The function C = 
of k = bw/v for harmonic motion. 
to calculate the forces in terms of h and a@ by standard pro- 
cedures and transform the results to w and @ coordinates. 


If 


C(k) is given in reference 1 as a function 
It is sometimes convenient 


= Z(a) + Z(h) 
M(a) + M(h) 


for harmonic h and a of frequency 


SN 


are the force and moment 
w, then Eqs. (5) through (8) are simply 


Z(w) = Z(h)/iw (5a) 
Z(0) = Za) — (u/iw)Z(h) (6a) 
M(w) = M(h)/iw (7a) 
M(@) = M(a) — (u/iw) M(h) (8a) 


Eqs. (5) and (7) could be obtained from Eqs. (1) and (2) by 
taking w = h and 6 = & = O, and Eqs. (6) and (8) could be 
obtained from Eqs. (1) and (2) by the condition A = w/u = 0 
fi.e.,0 = A = (h/u) + 0 which gives h = —u@andh = —u6 = 
This gives the physical interpretation of the wing moving 
Eqs. (5) and 


—va). 
in one of either the w or @ coordinate directions. 
(7) have been given in a number of previous works, but the @ 
expressions have been incorrect. The correct equations [Eqs. 
(6) and (8) above] reduce to those given by Glauert? for motion 
in which 6 = g = constant. 

In reference 2, Glauert also gives the effect of finite aspect 
ratio on the steady-state value. A comparison of Glauert’s 
value with the ordinary effect of aspect ratio on the lift curve 
slope, Crg, is shown in Fig. 2, and it is seen that a satisfactory 
approximation is obtained, at least for small values of k, by 
replacing 27 by Cr in the aerodynamic circulation terms de- 
termined by the velocity. The of Reissner* indicates 
that replacing 27 by Cz, will be satisfactory for the aerodynamic 
circulation terms only for k < 0.1, which range fortunately covers 
most airplane dynamics including the short period. At higher 
frequencies, the effect of aspect ratio is practically eliminated, 
but the error in the above approximation is not serious, since the 
aerodynamic circulation terms become negligible in comparison 
with the apparent mass terms depending on the accelerations. 

The apparent mass terms related to w and g for any plan form 
of small aspect ratio can be approximated by replacing pb? by 
p(2/3)(SC/E), which is the exact apparent mass factor for a 
finite flat elliptic plate in potential flow (see reference 4). The 
factor E depends only on the aspect ratio of the plate and is the 
complete elliptic integral of the second kind, whose numerical 
value is plotted in Fig. 2 and is the ratio of the elliptic plate 
Using these aspect ratio correc- 


work 


perimeter to twice the span. 
tions, Eq. (6) may be written in terms of wing area S and root 
chord C for rotation about the center of gravity located G 
distance from the leading edge: 


se, ‘ ° of . 

ne & _ ©) cigC()a so xP “ ( — ay (9) 
where Cz, is the static value of 0C,/da. 

In the same way, similar expressions corresponding to Eqs. 
(5), (7), and (8) may be written. 

It is interesting to note that the lift due to gq, the first term in 
Eq. (9), will produce a downwash even though A = w/u = 0. 
This occurs because the chordwise vortex distribution is changed 
to adjust effectively the wing camber to the streamline curvature 
produced by g, and its magnitude may be obtained as follows: 

It is reasonable to assume that the static value of de/dCz 
holds for nonsteady motion if only the circulatory part of the lift 
force is considered, since the circulation is the origin of the 
A change in the downwash angle e can be written as 


Z(0) 


downwash. 
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de = (de/dCy,) dCi 
and for a fixed velocity, u, 
dZ = dL = (pu?/2)S dCi 
Thus, 
de = (de/dCy)(2/pu2S) dZ 
Now, as discussed above, de/dC;, will be assumed to apply 
to the circulation part of the Z force developed in oscillatory 
motion, which, of course, depends on frequency. Let de 
and dZ be changes from the values at the midpoint of the oscilla- 
tion. 
Then, 
= (de/dCz)(2/pu?S)Z.(w) 
and 
€(0) = (de/dCz)(2/pu?2S)Z.(6) 
where the subscript c indicates the circulation part of Z |i.e., 
the part multiplied by the C(k) of reference 1]. 
This downwash occurs when the lift changes, and it will reach 
a point / distance behind the wing (e.g., the tail of a conventional 
airplane) after an interval Af = //u. For harmonic motion, 


the interval //u corresponds to a phase angle /w/u rad.; therefore 
the expressions above must be multiplied by e~*@e/"), 

At low frequencies (k < 0.1), the downwash at the tail iS not 
appreciably changed from the quasi-stationary value, and, ig 
particular, that due to q is negligible. Also, it is readily shopy 
that, for w < u/l, the phase lag factor e~*/¢/” approaches the 
downwash lag used in the conventional stability derivatives, 
With increasing frequency, these effects become of considerable 
importance up to approximately k’ = b’w/u > 27 (b' = chopg 
of tail). When the wave length of downwash variation is equal 
to, or less than, the tail chord, the net force on the tail approaches 


zero. 
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